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?
Thanks in advance.
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]$.
 A: 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
A: 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.
