function [p,L,px]=lagrange_interp1(x,f,varargin)
% [P,L,PX]=LAGRANGE_INTERP1(X,F,XX)
% DADAS UNAS ABSCISAS X=[X(1),...,X(N)], LOS CORRESPONDIENTES VALORES DE
% LA FUNCION F=[F(1),...,F(N)] Y UN CONJUNTO DE M PUNTOS XX=[XX(1),...,
% XX(M)]; CALCULA LOS POLINOMIOS DE LA BASE POLINOMICA DE LAGRANGE: L(1),
% ...,L(N); EL CORRESPONDIENTE POLINOMIO INTERPOLADOR GRADO N-1, P, Y LOS
% VALORES DE P EN LOS PUNTOS XX(1),...,XX(M).
%
% INPUT:
%     X=[X(1),...,X(N)]: CONJUNTO DE ABSCISAS PARA CALCULAR EL POLINOMIO
%                        INTERPOLADOR.
%     F=[F(1),...,F(N)]: VALORES DE LA FUNCION F EN X(1),...,X(N) RESPECT.
%  XX=[XX(1),...,XX(M)]: (OPCIONAL). CONJUNTO DE PUNTOS SOBRE EL EJE X EN
%                        LOS QUE SE EVALUARA EL POLINOMIO INTERPOLADOR.
% OUTPUT:
%     P=[P(1),...,P(N)]: COEFICIENTES DEL POLINOMIO INTERPOLADOR EN ORDEN
%                        DECRECIENTE DEL GRADO, QUE ES COMO DE ESCRIBEN
%                        LOS POLINOMIOS EN MATLAB/OCTAVE.
% L=[L(:,1),...,L(:,N)]: MATRIZ QUE CONTIENE, POR COLUMNAS, LOS COEFS. DE
%                        LOS POLINOMIOS QUE FORMAN LA BASE POLINOMICA DE
%                        LAGRANGE. ASI, L(:,1) CONTIENE LOS COEFS. DE L1,
%                        L(:,2) LOS DE L2,...,L(:,N) LOS DE LN.
%  PX=[PX(1),...,PX(M)]: VALORES DEL POLINOMIO INTERPOLADOR EN LOS PUNTOS
%                        XX(1),...,XX(M).

%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

% LAGRANGE_INTERP1
% COPYRIGHT (C) 2014 DEPT. MATEMATICA APLICADA I, UPC
% CREATED: 2014-02-20

n=length(x);
m=n-1;
L=zeros(n);
L(1:m,:)=[1:m]'*ones(1,n)+tril(ones(m,n));

for i=1:n
  p=poly( x( L(1:m,i) ) ); %p=(x-x(1))* * *(x-x(i-1))*(x-x(i+1))* * *(x-x(n))
  L(:,i)=p/polyval(p,x(i));%L_i=p(x)/p(x(i))
end
p=L*f'; p=p';

if length(varargin) > 0
   px=polyval(p,varargin{1});
end

return
end
