% function [y_int] = lagrange( x_int, x, y )
% Interpolate data by the Lagrange polynomial.
% Given points are (x,y), the polynomial is evaluated at points x_int.
function [y_int] = lagrange( x_int, x, y )
  
  n = length( x );
  np = length( x_int );
  
  lb = zeros( n, 1 );
  % Process all evaluation points.
  for p = 1:np
    % Create base.
    px = x_int( p );
    for i = 1:n
      jj = setdiff( 1:n, i );
      lb(i) = 1.0;
      for j = jj
	lb(i) = lb(i) * (px - x(j)) / (x(i) - x(j));
      end
    end
    % Interpolate.
    y_int( p ) = lb' * y(:);
  end