Intervals where the sign of a polynomial can be computed reliably

Question

This is a follow-up of a previous question.

Let $p$ be a polynomial with floating-point coefficients.

Is there a method for finding intervals where evaluating $p$ in floating-point arithmetic always gives the correct sign?

I want to find (ideally, maximal) disjoint interval $I_1, \dots, I_m$ such that for all floating-point numbers $a \in I_k$, the sign of $f(a)$ is always correct.

This is clearly related to isolating the zeros of $f$.

A standard example, where the sign of $f$ is wrong around a zero is $f(x)=x^6 - 6 x^5 + 15 x^4 - 20 x^3 + 15 x^2 - 6 x + 1$ near the multiple zero $x=1$. Note the $f(x)$ is the expansion of $(x-1)^6$.

Answer 1

Yes. You can compute a running error bound, i.e, a number $\mu$ such that the difference between the exact value of $y = p(x)$ and the computed value satisfies $\hat{y}$ satisfies $$|y - \hat{y}| \leq \mu u.$$ Here $u$ is the unit roundoff. You can trust the sign computed sign of $y$, when $|\hat{y}| > \mu u$.

Let $p(x) = \sum_{j=0}^n a_j x^j$, then Horner's method computes $$p_0 = a_n, \quad p_i = x p_{i-1} + a_{n-i}, \quad i = 1,2,\dotsc,n.$$ If $a_i$ and $x$ are machine numbers, then a running error bound can be computed as follows $$\mu_0 = 0, \quad z_j = p_{j-1} x, \quad p_j = z_j + a_{n-j}, \quad \mu_j = \mu_{j-1} |x| + |z_j| + |p_j|, \quad j=1,2,\dotsc, n.$$ The algorithm returns $y = p_n$ and $\mu = \mu_n$.

It is possible to reduce the cost of this algorithm from 5n flops to 4n flops, see Higham's book "Accuracy and Stability of Numerical Algorithms" for details.

Answer 2

I want to add that in addition to Carl Christian's suggestion of using a running error bound, you can also take the general relative error bound $$ \frac{|\hat p(x)-p(x)|}{|p(x)|} \leq \gamma_{2n}\,\mathrm{cond}(p,x),\qquad \gamma_{2n} \approx 2nu,\\ \mathrm{cond}(p,x) = \frac{\sum |a_i||x|^i}{|\sum a_i x^i|} = \frac{\tilde p(x)}{|p(x)|},$$ (see http://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf), and then the condition that the sign is right can be imposed be imposing $$ \gamma_{2n}\mathrm{cond}(p,x) \leq 1, \quad\text{or}\quad |p(x)|\gtrsim 2nu\tilde p(x)$$

Essentially all this is doing is identifying those $x$ for which evaluating the polynomial at $x$ is well-conditioned, so it's just a bound on the condition number.

For $(x-1)^6$, this would result in something like $|x-1| > c u^{1/6}$ where $c$ is a small constant $c$ and $u$ is the unit roundoff. For a simple root $\alpha$, it would be $|x-\alpha| > c u$.