function p = polint(xk, fk, x, alpxk, alpx)
% The function p = polint(xk, fk, x) computes the polynomial interpolant
% of the data (xk, fk). Two or more data points are assumed.
%
% Input (constant weight)
% xk: Vector of x-coordinates of data (assumed distinct).
% fk: Vector of y-coordinates of data.
% x: Vector of x-values where polynomial interpolant is to be evaluated.
%
% In the case of non-constant weights two additional input parameters are
% needed, as in p = polint(xk, fk, x, alpxk, alpx)
%
% Input (non-constant weight)
% xk: Vector of x-coordinates of data (assumed distinct).
% fk: Vector of y-coordinates of data.
% x: Vector of x-values where polynomial interpolant is to be evaluated.
% alpxk: Vector of weight values sampled at the points xk.
% alpx: Vector of weight values sampled at the points x.
%
% Output:
% p: Vector of interpolated values.
%
% The code implements the barycentric formula; see page 252 in
% P. Henrici, Essentials of Numerical Analysis, Wiley, 1982.
% (Note that if some fk > 1/eps, with eps the machine epsilon,
% the value of eps in the code may have to be reduced.)
% Note added May 2003: Except for certain nice node distributions
% polynomial interpolation of high-degree is an ill-conditioned
% problem. This code does not test for conditioning so use with
% care.
% J.A.C. Weideman, S.C. Reddy 1998.
if nargin == 3
alpx = 1; alpxk = 1;
elseif nargin > 3
fk = fk./alpxk;
end
x = x(:); % Make sure the data are column vectors
xk = xk(:); fk = fk(:);
alpxk = alpxk(:); alpx = alpx(:);
N = length(xk);
M = length(x);
L = logical(eye(N));
D = xk(:,ones(1,N))-xk(:,ones(1,N))'; % Compute the weights w(k)
D(L) = ones(N,1);
w = 1./prod(D)';
D = x(:,ones(1,N)) - xk(:,ones(1,M))'; % Compute quantities x-x(k)
D = 1./(D+eps*(D==0)); % and their reciprocals.
p = alpx.*(D*(w.*fk)./(D*w)); % Evaluate interpolant as
% matrix-vector products.