function [y, delta] = polyval(p,x,S,mu)
//y = polyval(p,x) returns the value of a polynomial of degree n evaluated at x. The input argument p is a vector of length n+1 whose elements are the coefficients in descending powers of the polynomial to be evaluated.
//x --can be a matrix or a vector. In either case, polyval evaluates p at each element of x.
//[y,delta] = polyval(p,x,S) uses the optional output structure S generated by polyfit to generate error estimates delta.
//delta --is an estimate of the standard deviation of the error in predicting a future observation at x by p(x).
//               If the coefficients in p are least squares estimates computed by polyfit, and the errors in the data input to poly//               if it are independent, normal, and have constant variance, then y±delta contains at least 50% of the
//                  predictions of future observations at x.

//y = polyval(p,x,[],mu) or [y,delta] = polyval(p,x,S,mu) use ˆx=(x−μ1)/μ2 in place of x. In this equation, μ1=mean(x)     //      and μ2=std(x). The centering and scaling parameters mu = [μ1,μ2] is optional

//EXAMPLES:
//p = [3 2 1];
//y=polyval(p,[5 7 9])
//EXPECTED OUTPUT:
//y=  86  162  262

if ~(isvector(p) | isempty(p))  // Check input is a vector
    error('polyval:InvalidP');
end

nc = length(p);
if isscalar(x) & (argn(2) < 3) & nc>0 & (abs(x)<%inf) & and(abs(p(:))<%inf)
    // Make it scream for scalar x.  Polynomial evaluation can be
    // implemented as a recursive digital filter.
    y = filter(1,[1 -real(x)],p);
    y = y(nc);
    return
end

siz_x = size(x);
if argn(2) == 4
   x = (x - mu(1))/mu(2);
end

// Use Horner's method for general case where X is an array.
y = zeros(size(x,1),size(x,2));
if nc>0, y(:) = p(1); end
for i=2:nc
    y = x .* y + p(i);
end

if argn(1) > 1
    if argn(2) < 3 | isempty(S)
        error('polyval:Requires S');
    end

    // 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('polyval:SizeS');
      end
      R = S(1:nc,1:nc);
      df = S(nc+1,1);
      normr = S(nc+2,1);
    end

    // Construct Vandermonde matrix for the new X.
    x = x(:);
    V(:,nc) = ones(length(x),1,class(x));
    for j = nc-1:-1:1
        V(:,j) = x.*V(:,j+1);
    end

    // S is a structure containing three elements: the triangular factor of
    // the Vandermonde matrix for the original X, the degrees of freedom,
    // and the norm of the residuals.
    E = V/R;
    e = sqrt(1+sum(E.*E,2));
    if df == 0
        warning('polyval:ZeroDOF');
        delta = %inf(size(e));
    else
        delta = normr/sqrt(df)*e;
    end
    delta = matrix(delta,siz_x);
end
endfunction