Jacques Gélinas on Sat, 09 Jun 2018 12:49:41 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

RE: Counting real roots of integer polynomials


The Budan-Fourier theorem is faster to use than Sturm's theorem.
It 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.
The bound is exact on every interval if and only if the polynomial has only real roots.

Here are one-liners showing the problems with multiple roots at the ends
(last two calls of "nsturm" in "ckBFS").

To avoid the large numbers appearing in Sturm sequences evaluations:
1. Apply "sturml" to see if the l.c. have the same signs (=>real roots only)
2. Use "nbudan" to locate the roots.

eva(v,a)=subst(v,variable(v),a); \\ one variable evaluation

sturml(p) = apply(pollead, sturmv(p));
addhelp(sturml,"sturml(p): leading coefficients of Sturm sequence for p");

nsturm(p,a,b) = my(stp=sturmv(p)); nvar(eva(stp,a)) - nvar(eva(stp,b));
addhelp(nsturm,"nsturm(p,a,b): number of distinct zeros of polynomial p in (a,b]");

nvar(V) = my(r,s); sum(k=1,#V,if(s,r=s);s=sign(V[k]);r&s&(r!=s));
addhelp(nvar,"nvar(V): number of variations of nonzero signs in real vector V");

sturmv(p,K) ={  \\ [a*x^2+b*x+c, 2*(a*x+b), (b^2-a*c)/a]
  my( m=poldegree(p)+1, stv=vector(m), k=2, stk=deriv(p) );
  stv[1]=p; while(stk&&(k<=(if(!K,m,K))), stv[k]=stk; k++; stk=-stv[k-2]%stv[k-1]);
  return(if(!K,stv,stv[1..K]));
}
addhelp(sturmv,"sturmv(p,{K}): first K (all by default) polynomials of Sturm sequence for p");

nbudan(p,a,b) = {
  local( BFv );
  if(p==0, error("nbudan: zero polynomial"));
  BFv = concat(p,vector(poldegree(p),k,p=p'));
  nvar(eva(BFv,a)) - nvar(eva(BFv,b));
}
addhelp(nbudan,"nbudan(p,a,b): Budan-Fourier root bound for polynomial p on (a,b]");

ckBFS() = { my( p = (x^2-1)*(x-2)^2, q = (x^2+1)*(x-2)^2 );
[ nvar(sturml(p)) == 0, nvar(sturml(q)) == 1, 
  nsturm(p,1,2) == 1, nsturm(p,0,2) == 2,
  nbudan(p,1,2) == 2, nbudan(p,0,2) == 3,
  nbudan(q,1,2) == 2, nbudan(q,2,2) == 0,
  nsturm(q,1,2) == 2,  \\ wrong !
  nsturm(q,2,3) == -1  \\ very bad !!
 ];}


Jacques Gélinas