function p = ipolevalbarycentric(x,y,t)
% Evaluation of interpolating polynomial based on barycentric 
% interpolation formula.
% x is a row vector of nodes, y a row vector of data, and t a row 
% vector of evaluation points
n = length(x); % number of interpolation nodes = degree of polynomial -1
N = length(t); % Number of evaluation points stored in t
% Precompute the weights
for k = 1:n
  lambda(k) = 1 / prod(x(k) - x([1:k-1,k+1:n])); end;
for i = 1:N
% Compute quotient of weighted sums 
  z = (t(i)-x); j = find(abs(z) <= eps*norm(z));  
  if (~isempty(j)), p(i) = y(j); % avoid division by zero
  else 
    mu = lambda./z; p(i) = sum(mu.*y)/sum(mu);
  end
end
