Budan-Fourier no-root test

[Jacques Gélinas <jacquesg00@hotmail.com>]

0. To prove inequalities, I need the sign of polynomials with integer coefficients on an interval with rational endpoints, while avoiding floating-point arithmetic. For example, the sign of 8*P^2 - 30*P + 15 is negative at P=Pi because it is negative on [3,22/7], and say I want these signs,  where u is the variable and P represents Pi (which is not algebraic) ? 
                            
psg((8*P^2-30*P+15)*u  + 3,P,8) == +1
psg(64*u^2 + (32*P^2 - 144*P + 210)*u + 96*P^2 - 240*P + 150 - 1200,P,oo) == +1
psg(64*u^2 + (32*P^2 - 144*P + 210)*u + 96*P^2 - 240*P + 150 - 1200,3,P) == Unk

1. The Budan-Fourier theorem bounds the number of roots with multiplicities of a real polynomial in a left-open interval by comparing the numbers of sign variations of the sequence of its derivatives evaluated at the ends of the interval. It generalizes the Descartes' rule of signs (on which its proof can be based {Conkwright:1943}, and has been extended to real analytic functions {Hurwitz:1912}. The bound is exact on every interval if and only if the polynomial has only real roots {Curtiss:1918}.

\\-------------------------------------------------------- begin cut
alias(dg,poldegree); 

nbudan(p,a,b) = {  local( BFv );
  if(p==0, error("nbudan: zero polynomial"));
  BFv = concat(p,vector(dg(p,u),k,p=deriv(p,u)));
  nvar(BFv,a) - nvar(BFv,b);
}

\\ 2. Variations of nonzero signs in real vector V.

nvar(V,a)=local(r,s); sum(k=1,#V, if(s,r=s); s=signa(V[k],a); if(r&&s,r!=s) );
signa(p,a)= signpi( if(a==oo, pollead(p,u), subst(p,u,a)) );

\\ 3. For the sign of coefficients involving P=Pi, recursion and substitution are used here.

pi=333/106; PI=355/113;
signpi(p) = if( dg(p,P)<=0, sign(simplify(p)), psg(subst(p,P,u),pi,PI) );

\\ 4.  The sign determination routine returns Unk if the  polynomial p possibly changes sign on (a,b).

psg(p,a=P,b=oo) = {  local(sga, sgb);
  sga = signa(p,a);
  if(dg(p,u)<=0||a==b, return(sga));
  sgb = signa(p,b);
  if(!sga||!sgb||sga!=sgb||nbudan(p,a,b), sga=Unk );  
  return(sga);
}
\\-------------------------------------------------------- end paste

5. Here are some problems with these routines.
a) The variable name u is fixed, which leads to  psg(subst(p,P,u),pi,PI)), not elegant.
b) Recursion should be avoided if the coefficients contain Pi . How ?
c) Also I would prefer to use more advanced intrinsic PARI/GP functions here.

6. Note that polsturm does not work with two variables (and is not efficient for a no-root test):
# polsturm((8*P^2-30*P+15)*u+3,3,8)
  *** polsturm: incorrect type in gsigne (t_POL).

Thanks for any suggestions for simplifications,

Jacques Gélinas, Ottawa