Most efficient way to compute a polynomial

Question

Polynomial: a0x^0 + a1x^1 +a2x^2 + a3x^3 + ... + anx^n

Array: array_a[] = {a0, a1, a2, a3 ... an};

I wrote a function to calculate this polynomial in Java:

public double cal(double x) {
    double y = 0.0;
    for (int index = array_a.length - 1; index >= 0; index--) {
        y = array_a[index] + y * x;
    }
    return y;
}

This seems 5 times faster than the loop y += array_a[index] * Math.Pow(x, index);

But I wondering if there is a better way to compute this polynomial?

** For anyone thinks it's a different calculation: I did test the function above. It does the same thing with y += array_a[index] * Math.Pow(x, index); and they compute the same result.

Thanks.

Answer 1

This is Horner's method. If you only want to calculate it once per polynomial, this is the most efficient algorithm:

… Horner's method requires only n additions and n multiplications, and its storage requirements are only n times the number of bits of x. …

Horner's method is optimal, in the sense that any algorithm to evaluate an arbitrary polynomial must use at least as many operations. Alexander Ostrowski proved in 1954 that the number of additions required is minimal. Victor Pan proved in 1966 that the number of multiplications is minimal.

If you need to evaluate the polynomial extremely many times and the degree is very high, then there are methods to transform the representation of the polynomial (preconditioning) so that the number of multiplication is reduced to ⌊n/2⌋ + 2. This seems not very practical though, at least I've never seen this in the wild. I've found an online paper that describes some of the algorithms if you are interested.

Also mentioned in the paper, due to the CPU architecture it might be more efficient if you evaluating even and odd terms separately so they can be placed in parallel pipelines:

public double cal(double x) {
    double x2 = x * x;
    double y_odd = 0.0, y_even = 0.0;
    int index = array_a.length - 1;
    if (index % 2 == 0) {
        y_even = array_a[index];
        index -= 1;
    }
    for (; index >= 0; index -= 2) {
        y_odd = array_a[index] + y_odd * x2;
        y_even = array_a[index-1] + y_even * x2;
    }
    return y_even + y_odd * x;
}

The JIT/compiler might be able to do this conversion for you or even use SIMD to make it very fast automagically. Anyway, for this kind of micro-optimization, always profile before committing to a final solution.