Algorithm to check if a polynomial is positive in a given interval

I'm writing a program where I have a 3-dimensional polynomial (3 variables) of which I have to check if it is positive in a given product of 3 intervals (a volume).

I found a paper which solved this problem but I didn't really understand it, nor am I able to extract an algorithm out of it: http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1084864

So can this approach be formulated as an algorithm for a computer?

Edit:

To make things more specific.

My polynomial is always:

z^3*(z*(6*z-15)+10)(-y^3(y*(6*y-15)+10)(x^3(x*(6*x-15)+10)*(g3z*z-g2z*z+g3y*(y-1)-g2y*(y-1)-g2x*x+g3x*(x-1))-x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))+g2z*z-g0z*z-g0y*y+g2y*(y-1)+g2x*x-g0x*x)-x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))-g0z*z+g4z*(z-1)+y^3*(y*(6*y-15)+10)(g6z(z-1)-g4z*(z-1)+x^3*(x*(6*x-15)+10)(g7z(z-1)-g6z*(z-1)+g7y*(y-1)-g6y*(y-1)-g6x*x+g7x*(x-1))-x^3*(x*(6*x-15)+10)(g5z(z-1)-g4z*(z-1)+g5y*y-g4y*y-g4x*x+g5x*(x-1))-g4y*y+g6y*(y-1)+g6x*x-g4x*x)+x^3*(x*(6*x-15)+10)(g5z(z-1)-g4z*(z-1)+g5y*y-g4y*y-g4x*x+g5x*(x-1))+g4y*y-g0y*y+g4x*x-g0x*x)+y^3*(y*(6*y-15)+10)(x^3(x*(6*x-15)+10)*(g3z*z-g2z*z+g3y*(y-1)-g2y*(y-1)-g2x*x+g3x*(x-1))-x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))+g2z*z-g0z*z-g0y*y+g2y*(y-1)+g2x*x-g0x*x)+x^3*(x*(6*x-15)+10)*(g1z*z-g0z*z+g1y*y-g0y*y-g0x*x+g1x*(x-1))+g0z*z+g0y*y+g0x*x

where all variables starting with $g$ are constant. The given intervals for $x,y,z$ are limited to $[0,1]$.

-
Do you really mean a $3$-dimensional polynomial, or a polynomial of degree $3$? If by a $3$-dimensional polynomial you mean a polynomial containing $3$ variables, then it doesn't make sense to talk of one given interval. That article doesn't seem to be freely available to non-members. – joriki Jul 10 '11 at 22:34
You might want to look up Sturm sequences. – Yuval Filmus Jul 10 '11 at 22:40
I'm talking about a 3-dimensional polynomial (means 3 variables) and a 3-dimensional interval which is a volume (en.wikipedia.org/wiki/…). The degree is about 6. – man Jul 10 '11 at 22:44
@man: "3-dimensional interval" is not a standard use of the word "interval." If you mean a product of intervals, say so. – Qiaochu Yuan Jul 10 '11 at 22:50
Interesting -- "n-dimensional interval" does turn up about 50,000 Google hits -- I'd never seen that expression before. – joriki Jul 10 '11 at 23:00

The paper you cite contains the following theorem:

$$P(x) > 0, \ \ \forall x \in D^n \subseteq R^n$$

where $P(x)$ is a real polynomial, iff:

1) All the $n-1$-dimensional polynomials arrived at by substituting the end points of the intervals of $x_i$ $(i=1,2,\cdots,n)$ in $P(x)$ (each one at a time), are positive in $D^{n-1}$,

2) The set of $n$ equations in $n$ real variables \eqalign{P(x) &= 0 \cr \frac{\partial P}{\partial x_i} &= 0, \ \ \forall i \in \{1,2,\cdots,n-1\}\cr} has no solution in $D^n$.

I would program this in Maple (for the case of finite intervals) as follows:

ispositive:= proc(P::polynom, R::list(name=range))
local x,D,n,i,S,s,Q;
x:= map(lhs,R);
D:= map(convert,map(rhs,R),list);
if D = [] then return is(P > 0) end if;
n:= nops(D);
for i from 1 to n do
if not (ispositive(subs(x[i]=D[i][1], P),subsop(i=NULL,R))
and ispositive(subs(x[i]=D[i][2], P), subsop(i=NULL,R))) then return false
end if
end do;
S:= RootFinding[Isolate]({P, seq(diff(P,x[i]),i=1..n-1)},x,output=interval);
for s in S do
Q:= map(rhs,s);
if or(seq(Q[i][1] > D[i][2],i=1..n),seq(Q[i][2] <D[i][1],i=1..n)) then next end if;
if and(seq(Q[i][1] >= D[i][1],i=1..n),seq(Q[i][2]<=D[i][2],i=1..n))
then return false
else error "Roots returned with too little precision, try increasing Digits"
end if;
end do;
true
end proc;


For example,

ispositive(16*x^2 + 16*y^2 - 64*x - 16*y + 68,[x=1..10, y=-1..1]);

false

-
Wow, thank you. But what still confuses me is that you are using RootFinding[Isolate] which only works, if there is a finite number of solutions. – man Jul 11 '11 at 20:10
Yes, so this may not work in some "degenerate" cases. For example, ispositive((x^2+y^2-1)^2,[x=-2..2,y=-2..2]); – Robert Israel Jul 12 '11 at 0:31

I solved this by converting the polynomial to a bezier volume with regular control points and splitting it to lie in my asked intervals. Then I can use the fact that the curve lies completely inside the convex hull of its control points.

-
This is a very good approach to the problem. In general you can use Interval Arithmetic (along with some subdivision) to solve this sort of thing, but the right form for the polynomial can make a huge difference in the results. The convex-hull property of Bezier surfaces that you mention is a perfect example of this. – Steven Stadnicki Aug 24 '11 at 17:08