function [y, delta] = polyval(p,x,S)
%POLYVAL Evaluate polynomial.
%   Y = POLYVAL(P,X), when P is a vector of length N+1 whose elements
%   are the coefficients of a polynomial, is the value of the
%   polynomial evaluated at X.
%
%       Y = P(1)*X^N + P(2)*X^(N-1) + ... + P(N)*X + P(N+1)
%
%   If X is a matrix or vector, the polynomial is evaluated at all
%   points in X.  See also POLYVALM for evaluation in a matrix sense.
%
%   [Y,DELTA] = POLYVAL(P,X,S) uses the optional output structure
%   generated by POLYFIT to generate error estimates, Y +/- delta.
%   If the errors in the data input to POLYFIT are independent normal
%   with constant variance, Y +/- DELTA contains at least 50% of the
%   predictions.
%
%   See also POLYFIT, POLYVALM.

%   Polynomial evaluation p(x) using Horner's method

%   J.N. Little 4-26-86
%   Revised 3-9-88 JNL, 5-5-93 BAJ
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.10 $  $Date: 1997/11/21 23:40:58 $

[m,n] = size(x);
nc = length(p);
infs = find(isinf(x));
sinfs = (sign(x(infs)).^(length(p))).*x(infs); 

if ((m+n) == 2) & (nargin < 3)
    % Make it scream for scalar x.  Polynomial evaluation can be
    % implemented as a recursive digital filter.
    y = filter(1,[1 -x],p);
    y = y(nc);
    return
end

% Do general case where X is an array
y = zeros(m,n);
for i=1:nc
    y = x .* y + p(i);
end

if nargin > 2 & nargout > 1
    x = x(:);

    % Extract parameters from S
    if isstruct(S),  % Use output structure from polyfit.
      R = S.R;
      df = S.df;
      normr = S.normr;
    else             % Use output matrix from previous versions of polyfit.
      [ms,ns] = size(S);
      if (ms ~= ns+2) | (nc ~= ns)
          error('S matrix must be n+2-by-n where n = length(p)')
      end
      R = S(1:nc,1:nc);
      df = S(nc+1,1);
      normr = S(nc+2,1);
    end

    % Construct Vandermonde matrix.
    V(:,nc) = ones(length(x),1);
    for j = nc-1:-1:1
        V(:,j) = x.*V(:,j+1);
    end

    % S is a structure containing three elements: the Cholesky factor of the
    % Vandermonde matrix, the degrees of freedom and the norm of the residuals.
    E = V/R;
    if nc == 1
        e = sqrt(1+(E.*E));
    else
        e = sqrt(1+sum((E.*E)')');
    end
    if df == 0
        disp('Warning: zero degrees of freedom implies infinite error bounds.')
        delta = repmat(Inf,size(e));
    else
        delta = normr/sqrt(df)*e;
    end
    delta = reshape(delta,m,n);
end

y(infs) = sinfs;
