function [y,delta] = polyval(p,x,S,mu)
%MATLAB Code Generation Library Function

%   Limitations:
%   Matrix S input from previous versions of POLYFIT is not supported.

%   Copyright 1984-2014 The MathWorks, Inc.
%#codegen

eml_invariant(nargin >= 2, ...
    'MATLAB:minrhs');
eml_invariant(isvector(p), ...
    'MATLAB:polyval:InvalidP');
eml_invariant(isa(p,'float'), ...
    'Coder:toolbox:unsupportedClass','polyval',class(p));
eml_invariant(isa(x,'float'), ...
    'Coder:toolbox:unsupportedClass','polyval',class(x));
if nargin == 4
    eml_invariant(isa(mu,'float') && numel(mu) >= 2, ...
        'Coder:toolbox:polyval_invalidMu');
    if isa(x,'double') && isa(mu,'single')
        % Must convert x to single.
        [y,delta] = polyval(p,single(x),S,mu);
        return
    end
    if isreal(x) && ~isreal(mu)
        % Must make x complex.
        [y,delta] = polyval(p,complex(x),S,mu);
        return
    end
    x = coder.internal.div(x-mu(1),mu(2));
end
y = coder.nullcopy(eml_expand(coder.internal.scalarEg(x,p),size(x)));
nc = numel(p);
if ~(isempty(y) || isempty(p))
    % Use Horner's method.
    y(:) = p(1);
    for k = 2:nc
        y = (x .* y) + p(k);
    end
end
if nargout > 1
    eml_invariant(nargin >= 3 && ~isempty(S), ...
        'Coder:toolbox:polyval_6', ...
        'IfNotConst','Fail');
    eml_invariant(isstruct(S), ...
        'Coder:toolbox:polyval_7', ...
        'IfNotConst','Fail');
    % 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.
    % What follows is a streamlined algorithm for evaluating:
    %     E = V/R;
    %     e = sqrt(1+sum(E.*E,2));
    %     delta = normr/sqrt(df)*e;
    delta = coder.nullcopy(eml_expand(coder.internal.scalarEg(x,S.R),size(x)));
    ZERO = coder.internal.scalarEg(x,S.R);
    normrddf = coder.internal.div(S.normr,sqrt(S.df));
    if S.df == 0
        if ~coder.setwarningflag
            coder.internal.warning('Coder:MATLAB:polyval_ZeroDOF');
        end
        delta(:) = inf;
    elseif nc == 0
        delta(:) = normrddf;
    else
        w = coder.nullcopy(eml_expand(coder.internal.scalarEg(x,S.R),[nc,1])); % Work vector.
        ncm1 = nc - 1;
        nx = numel(x);
        for k = 1:nx
            w(nc) = 1;
            for j = ncm1:-1:1
                w(j) = w(j+1)*x(k);
            end
            s2 = ZERO;
            for j = 1:nc
                wj = coder.internal.div(w(j),S.R(j,j));
                s2 = s2 + wj*wj;
                for i = j+1:nc
                    w(i) = w(i) - wj*S.R(j,i);
                end
            end
            delta(k) = normrddf * sqrt(1 + s2);
        end
    end
end
