function [yp] = polyinterp(xd,yd,xp)
% computes the values of P(x) at the points xp
% where P is the Lagrange interpolating polynomial
% for the data set (xd,yd)

n=length(xd);
yp=zeros(size(xp));

for i=1:n % summation loop
   Lprod=ones(size(xp));
   for k=[1:i-1 i+1:n] % product loop
       % product by repeated multiplication
       Lprod=Lprod.*(xp-xd(k))./(xd(i)-xd(k))
   end
   yp=yp+yd(i)*Lprod; % summation by repeated addition 
end

end

