polynomial evaluation accuracy, multiplication versus division

Question

let us say I have have polynomial in x, divided by a power of x:

p = (a + x(b + x(c + ..)))/(x**n)

efficiency aside, which would be more accurate computation numerically, the above or using division:

p = (((a/x + b)/x + c)/x + ...)

thanks

Answer 1

In theory, there shouldn't be any difference - if the values are calculated accurately with 'infinite' precision.

Kernighan and Plauger state in their antique but excellent book 'Elements of Programming Style', that:

A wise programmer once said, "Floating point numbers are like little piles of sand; every time you move one, you lose a little sand and gain a little dirt".

The division has slightly fewer operations overall, which means there is slightly less opportunity to lose sand and gain dirt.

A detailed analysis would probably require a look at the coefficients (a, b, c, etc) as well perhaps as the value of x - what works when x is huge may not work well when x is close to zero, nor vice versa.

Answer 2

I think the difference is minimal, unless there is a chance that x**n overflows or underflows, in which case you should use the second expression.

The two expressions differ in two places:

1. The evaluation order is reversed (..., c, b, a) for the first expression and (a, b, c, ...) for the second expression. Which one is best depends on the value of the coefficients.
2. The first expression has the .../x**n at the end. As Jonathan explains, for that reason it may be expected that the second expression is more accurate, because it has fewer operations. However, I think that the .../x**n causes only a minimal loss of accuracy (compared to other places where you lose accuracy), unless the x**n overflows or underflows.

Answer 3

The answers provided are unfortunately wrong.

The second equation p = (((a/x + b)/x + c)/x + ...) is only marginal worse for accuracy and much, much worse for speed.

Why ? The relative errors for multiplication has only the main linear term and a small quadratic term. Division in contrast introduces higher, but very small terms (cubic, quartic):

e = relative error, assumed constant for both terms

a*b = a(1+e)*b(1+e) = a*b (1+2e+e^2) // multiplication

a/b = a(1+e)/b(1+e) = a/b (1+e)(1+e+e^2+e^3+...geometric series) // division

So division is always a bit worse than multiplication. For speed considerations: Divisions are always slower than multiplications, the normal factor can vary from 3x - 10x. So nested divisions are much slower than nested multiplications if you don't compute the last factor x^n not by pow(), but by nested multiplication.

The x^n can be easily computed by a loop multiplying the result double power = x; for (n-1) power *= x;

If you use pow(), be aware that it is mostly conveniently computed by exponential and logarithm, taking much more time than necessary (100x).

Are you aware that while the error between double and exact result remains small, polynominal results are very sensitive to changes of x for higher n's ?! So if you use higher n's be aware that your answers may be totally off the mark because small errors in x are astronomically amplified.