6
$\begingroup$

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

Is there a method for accurately evaluating the sign of $p(a)$ in floating-point arithmetic?

I don't care about the actual value $f(a)$, only about its sign. But I want the exact sign.

If it helps, we can assume that $ 0 \le a \le 1$.

Improve this question
$\endgroup$
8
  • $\begingroup$ I don't see how Horner can be beaten here. $\endgroup$
    – J. M.
    Sep 1, 2017 at 15:06
  • $\begingroup$ @J.M., how does Horner find the sign? $\endgroup$
    – lhf
    Sep 1, 2017 at 15:07
  • 1
    $\begingroup$ You can use compensated Horner instead of normal Horner. And there are (robust) algorithms for computing the determinant (or its sign) of a companion matrix of the polynomial. $\endgroup$
    – gammatester
    Sep 1, 2017 at 18:53
  • 2
    $\begingroup$ Note that "Is $p(a)$ exactly zero?" is a sub-problem of your problem, so computing the exact sign is probably quite hard: it seems like you'd need at least $\mathrm{degree}\times(\text{mantissa bits})$ bits of precision to precisely distinguish the value of $p(a)$ from zero. $\endgroup$
    – Kirill
    Sep 1, 2017 at 20:04
  • 1
    $\begingroup$ If you really want the true sign at any cost, you could just convert all the floating-point numbers to exact rational numbers, and evaluate it with exact rational arithmetic. $\endgroup$
    – Kirill
    Sep 2, 2017 at 1:31

3 Answers 3

Reset to default
3
$\begingroup$

Compensated Horner method (http://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf) has an error bound of the form $$ |\mathrm{comphorner}(p, x) - p(x)| \leq u|p(x)| + \gamma_{2n}^2\tilde p(x), \qquad \tilde p(x) = \sum_i |a_i||x|^i, \quad \gamma_n = \frac{n u}{1-n u}, $$ where $u$ is the unit roundoff, so it will give the correct sign so long as $|p(x)| \geq \gamma_{2n}^2\tilde p(x)$. This will hold everywhere except very near the roots. Since you can scale your polynomial so that $|a_i|\leq 1$ also holds in addition to $|x|\leq 1$, the region of failure around a root is about the size $u^{2/m}$, where $m$ is the multiplicity of the root. For single roots this means failure is only possible if the given argument is the root exactly, and for double roots it must be with a few machine epsilons.

I don't think there's a way to do this that works in general without relying on extra precision, it's a very ill-conditioned problem near the roots.

Improve this answer
$\endgroup$
2
$\begingroup$

You could use interval arithmetic for this purpose.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Interval arithmetic will give an interval $[L,U]$ containing $f(a)$. If $L>0$ or $U<0$, then we'll know the sign of $f(a)$. But what if $L \le 0 \le U$? $\endgroup$
    – lhf
    Sep 1, 2017 at 16:50
  • 1
    $\begingroup$ Then it's inconclusive but at least in the other two cases you'd get a guarantee on the sign that's as good as a rigorous proof. $\endgroup$
    – Juan M. Bello-Rivas
    Sep 1, 2017 at 18:10
1
$\begingroup$

You can use expansion arithmetics [1] combined with arithmetic filters [2] to improve performance. The idea is to first do the computation with standard floating point arithmetics together with bounds indicating whether the sign is accurate, and only use the (more expensive) expansion arithmetic when the sign could not be evaluated. See also my survey article [3] together with the companion implementation in my programming library geogram [4]. It has a standalone implementation as a .h/.cpp pair of files (use multiprecision_PSM from the download section).

[1] https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf

[2] https://hal.inria.fr/inria-00344297/fr/

[3] https://hal.inria.fr/hal-01225202

[4] http://alice.loria.fr/software/geogram/doc/html/index.html

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ I already specified "my article" and "my project" when referring to my stuff ([3], [4]). Please tell me if it is not explicit enough. I have no relation with [1] and [2]. $\endgroup$
    – BrunoLevy
    Sep 2, 2017 at 13:31
  • $\begingroup$ I did not notice it. Sorry about that, I will remove my comment. $\endgroup$
    – nicoguaro ♦
    Sep 3, 2017 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.