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$.

$\endgroup$
  • $\begingroup$ I don't see how Horner can be beaten here. $\endgroup$ – J. M. Sep 1 '17 at 15:06
  • $\begingroup$ @J.M., how does Horner find the sign? $\endgroup$ – lhf Sep 1 '17 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 '17 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 '17 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 '17 at 1:31
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.

$\endgroup$
2
$\begingroup$

You could use interval arithmetic for this purpose.

$\endgroup$
  • 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 '17 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 '17 at 18:10
0
$\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

$\endgroup$
  • 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 '17 at 13:31
  • $\begingroup$ I did not notice it. Sorry about that, I will remove my comment. $\endgroup$ – nicoguaro Sep 3 '17 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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