function y = Polyvec(c,A,p)
% c is a column m-vector
% A is an n-by-n matrix with lower bandwidth p, assume pm<n
% y is the first column of c(1)I + c(2)A + ... + c(m)A^{m-1}

[n,n] = size(A);
d = length(c);

% A 2-point solution that treats A as a full matrix...
y = zeros(n,1); y(1) = c(d);
for k=d-1:-1:1
y = A*y;
y(1) = y(1) + c(k);
end

% A 4 point solution that exploits the fact that at the start of the jth pass
% through the loop, the only nonzero part of y is y(1:1+(j-1)*p)...

y = zeros(n,1); y(1) = c(d);
j = 1;
for k=d-1:-1:1
y(1:1+j*p) = A(1:1+j*p,1:1+(j-1)*p)*y(1:1+(j-1)*p);
y(1) = y(1) + c(k);
j=j+1;
end

% A full-credit solution that is the same as the preceding, except that it
% exploits the bandedness of A(1:1+j*p,1:1+(j-1)*p)...

y = zeros(n,1); y(1) = c(d);
j = 1;
for k=d-1:-1:1
%y(1:1+j*p) = A(1:1+j*p,1:1+(j-1)*p)*y(1:1+(j-1)*p);
w = zeros(1+j*p,1);
for i=1:1+(j-1)*p
w(1:p+i) = w(1:p+i) + A(1:p+i,i)*y(i);
end
y(1:1+j*p) = w;
y(1) = y(1) + c(k);
j=j+1;
end