# A Brief History of polyval

When I first started working at MathWorks, there were already a few functions in MATLAB for working with polynomials. One
of these, `polyval`, for polynomial evaluation, used one algorithm, based on Horner's rule, known for its efficiency, especially suitable for embedded controller environments. After a customer wondered why the function
took so long to evaluate a given polynomial for a scalar input value `x`, we updated `polyval` to use a faster algorithm for the scalar case, and that algorithm uses the function `filter`. I thought I'd discuss a few more details here.

Note: since I wrote this post (but before I published it) there was a discussion on the newsgroup about this same topic. Be sure to read the reply from `John D'Errico`.

### Contents

### Horner's Rule

Let's first look at the code for evaluating a polynomial for multiple values of input `x`. Here's the relevant portion of `polyval`, `nc` is the length of the vector representing the polynomial (one more than the polynomial degree).

dbtype polyval 62:67

62 % Use Horner's method for general case where X is an array. 63 y = zeros(siz_x, superiorfloat(x,p)); 64 if nc>0, y(:) = p(1); end 65 for i=2:nc 66 y = x .* y + p(i); 67 end

Here's what a general polynomial looks like:

This equation can be rewritten like this,

and is the basis for evaluating the polynomial for an array `x`. Notice that the formula basically starts with the coefficient of the highest degree (named ) multiplied by `x`, then adds the next coefficient, multiplies that result by `x` again, adds the next coefficient, and so on, until reaching the final coefficient (the constant term, for ). In this formulation, we never need to explicitly create all the powers of `x` to evaluate the polynomial.

### The Scalar Case, Using filter

Here's the code for the scalar evalation case using the function `filter`.

dbtype polyval 52:53

52 y = filter(1,[1 -x],p); 53 y = y(nc);

It may not, at first glance, appear relevant to polynomial evaluation. However, taking a close look at the mathematics reveals
that it indeed can be useful. First, let's examine the `filter` formula. In this version, the coefficients for `a` have been normalized so `a(1)` = 1 (which is true for the inputs to `filter` you've just seen).

Now, let's try a small example. Suppose we want to evaluate the same third degree polynomial I used in another set of posts recently.

p = [1 0 1 -1];

which represents

When you call `filter` like this `y = filter(b,a,x)`, the results are

Looking again at the `polyval` code, use these substitutions for the inputs to `filter`: `b = 1`, `a = [1 -x]`, and `x = p`. We evaluate the filter for the values `p`, the polynomial coefficients! Let's do a little math to see why `filter` works. Following the formula for the output of the function `filter`, and assuming `p(0) = 0`, we can find the first output.

And here are the rest, going up to the final coefficient in our example.

It sure looks like the final answer here is the polynomial `p` evaluated at `x`.

### How Close are the Algorithms?

Numerically, the two algorithms end up being essentially the same.

x = 0:0.001:1; resultArray = polyval(p,x); resultScalar = zeros(size(x)); for index = 1:numel(x) resultScalar(index) = polyval(p,x(index)); end allequal = isequal(resultArray,resultScalar)

allequal = 1

Identical!

### Timing

I mentioned we made this change for speed. This particular code is tricky to check timing since, especially if the polynomial degree is small, the number of mathematic operations can be swamped by the timing mechanism. When computers were slower, and before we increased performance of for-loops in MATLAB, measuring the performance differences between these code alternatives was much easier than it is today.

### Multiple Algorithms

Do you try to write your code as one size fits all, or do you create code in which you switch algorithms or implementations,
such as we've seen in `polyval`, to account for different conditions such as array size? I'd love to hear more about your experiences with these ideas here.

## 评论

要发表评论，请点击 此处 登录到您的 MathWorks 帐户或创建一个新帐户。