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.

Partridgeberry

He uses that the polynomial is equal to a0 + x * (a1 + x*(a2 +... )) – Gooey

@ErwinBolwidt Yes I did. I know it looks like a different calculation. But it does the same thing. You can check. – Partridgeberry

Is there a better way? It seems like the best way already. Minimal number of computations. Why are you asking? What part don't you like? – Gentoo

@Gentoo Yes, I think so. This is the best way but my professor still asked her students to optimize this. So I'm trying to find a better way. – Partridgeberry

Maybe he wants this double y = 0; double x_n = 1; for (i = 0; i < length; i++) { y += array[i] * x_n; x_n *= x; } - performance wise it's probably equivalent to your version. – Nernst

@Nernst Performance-wise, your suggestion is doing two multiply operations per loop, while OP is only doing one. Not knowing how JIT may affect it, your suggestion would appear to be slower. – Gentoo

@Gentoo you are right – Nernst

No matter how you fiddle with it, you have to perform n add operations and n+1 multiply operations. That cannot be reduced. By using the a0 + x * (a1 + x*(a2 +... )) formula, you have successfully eliminated the power operations. Well, technically, multiplying by x^0 can be eliminated, so only n multiply operations. Your code is currently doing n+1 adds and multiply's, so you could twiddle it slightly. Note that n = array_a.length - 1. – Gentoo

All you can really do is start with double y = array_a[array_a.length -1] and then only run the loop from array_a.length-2. But that's such a marginal gain that it's hardly worth the slightly reduced readability. – Teeny

By the way, for practical purposes I'd probably store the array of coefficients the other way around, with the 0th element being a0 and so on. – Teeny

@Teeny it is already stored as you said so a0 at the index 0 ... etc. @ Louis, do you have the roots of the polynomial ? or it is just a general case, what about factoring it ? is that possible? – Harwill

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.

Amphibrach