function [M]=gera_base_discr(x,m)
// [M]=gera_base_discr(x,m)
// Gera uma base ortogonal de polinomios no sentido  
// do produto interno <f,g> = sum_{x(i)} f(x(i))*g(x(i))
// A entrada é o vetor de pontos x, por ex. x = [1 2 3 4 5]
// e a saida é uma matriz cujas linhas são os coefs. 
// dos pols ortogonais, por exemplo se 
// M = [ 1 0 0 0 ; .5 1 0 0 ; 0 1 1 0 ; 0 1/3 0 1] 
// então a base é phi_0 = 1, phi_1 = .5 + x, 
// phi_2 = x + x^2,   phi_3 = x/3 + x^3.
// São gerados phi_0 até phi_m sendo m entrada da função.
//
// Eduardo Fontoura 03/Maio/2016

n = size(x,2);

// matriz com potencias de x
P=[];
for i = 0:m,
P = [ P; x^i]
end

//phi_0:
M = zeros(m+1,m+1);
M(1,1)=1;

//phi_1:
M(2,2)=1; 
paux1 = conv([0 1],M(1,:));  // multiplicando o polinomio x pelo phi_0
paux1 = paux1(1:m+1); // tirando elementos desnecessários (e nulos)
av_paux1 = paux1*P;  // polinomio x*phi_0 avaliado nos pontos x
av_paux2 = M(1,:)*P;  // polinomio phi_0 avaliado nos pontos x
M(2,1)= -(av_paux1 * av_paux2') / (av_paux2 * av_paux2');

//demais phi_k,  k=2,...,m:
for k = 2:m
  paux1 = conv([0 1],M(k,:));  // x * phi_{k-1}
  paux1 = paux1(1:m+1);   // tirando elementos desnecessários (e nulos)
  av_paux1 = paux1*P;  // x*phi_{k-1} avaliado nos pontos x
  av_paux2 = M(k,:)*P;  // polinomio phi_{k-1} avaliado nos pontos x
  B= (av_paux1 * av_paux2') / (av_paux2 * av_paux2');

  av_paux1 = paux1*P;  // x*phi_{k-1} avaliado nos pontos x
  av_paux2 = M(k-1,:)*P;  // polinomio phi_{k-2} avaliado nos pontos x
  C= (av_paux1 * av_paux2') / (av_paux2 * av_paux2');

  Maux = conv([-B 1],M(k,:))
  Maux = Maux(1:m+1); // tirando elementos desnecessários (e nulos)
  M(k+1,:) = Maux - C*M(k-1,:)
end
endfunction


function [z,H,b]=minqd_base_geral(x,y,base)
// z=minqd(x,y,base)
// realiza aproximacao por minimos quadrados 
// discretos. P_n(x)=z_1 phi_1+...+z_{n+1} phi_{n+1} 
// eh o polinomio aproximador para os pontos dados
// e cada phi eh dado via a matriz base da seguinte 
// forma. Por exemplo, considerando um pol P_4, 
// phi_3 = 1+2x+3x^2 eh representado pelo 
// vetor [ 1 2 3 0 0] na 3a linha da base; 
// phi_1 = 1 => [ 1 0 0 0 0 ] na 1a linha, 
// phi_5 = 2x -x^4 => [ 0 2 0 0 -1] a 5a linha.
//
// Eduardo F Costa 02/05/2016
//
n = size(base,1)-1;
m = size(x,2);
err = 1e-6; // erro para integracao

// matriz com potencias de x
aux=[];
for i = 0:n,
aux = [ aux x'^i]
end

// phis avaliados nos pontos x e empilhados,  
// um ao lado do outro:
aux = (base*aux')'

// matriz do sistema de equações normais:
H = aux'*aux 
b=(aux')*y'  

z=H\b
endfunction


function [p,M,z] = minqdortog(x,y,n)
// z = [z_0 ... z_n] tal que P_n(x)=z_0 phi0+
// z_1 phi1+...+z_n phin eh o polinomio
// minimizador de erro quadrado discreto
// nos pontos x=[x_0,...,x_n], y=[y_0,...,y_n]
// p sao os coeficientes na base convencional, 
// tais que P_n(x)=p_0+...+p_n z^n
//
// Eduardo F Costa 02/05/2016
//

// DEFININDO A BASE (PHIs)
deff('y=f0(x)','y=1*x^0');
deff('y=fsum(x)','y=x');
B(1)=sum(x.*f0(x).*f0(x))/sum(f0(x).*f0(x));
deff('y=f1(x)','y=x-B(1)');

for k=2:n
disp(k)
deff('y=fnum(x)','y=((f'+string(k-1)+'(x)).^2).*x');
deff('y=fden(x)','y=((f'+string(k-1)+'(x)).^2)');
B(k)=sum(fnum(x))/sum(fden(x));
deff('y=fnum(x)','y=(f'+string(k-1)+'(x)).*(f'+string(k-2)+'(x)).*x');
deff('y=fden(x)','y=(f'+string(k-2)+'(x)).^2');
C(k)=sum(fnum(x))/sum(fden(x));
// definindo o fi_k
B_=B(k); C_ = C(k);
deff('y=f'+string(k)+'(x)','y=(x-'+string(B_)+').*(f'+string(k-1)+'(x))-'+string(C_)+'.*(f'+string(k-2)+'(x))');

end


// CALCULANDO A MATRIZ DE HILBERT
for i=0:n
for j=0:n
 deff('y=fmult(x)',strcat(['y=(f',string(i),'(x)).*(f',string(j),'(x))']));
 M(i+1,j+1)=sum(fmult(x));
end
end
disp('Veja como eh a estrutura da matriz de Hilbert:')
disp(M)


// CALCULANDO OS COEFICIENTES
for i=0:n
 deff('y=fmult(x)',strcat(['y=y.*(f',string(i),'(x))']));
 z(i+1)=sum(fmult(x))/M(i+1,i+1);  
end


deff('y=P0(x)','y=z(1).*x.^0');
for i=1:n-1
deff('y=P'+string(i)+'(x)','y=(P'+string(i-1)+'(x))+(z(i+1).*(f'+string(i)+'(x)))');
end
deff('y=P(x)','y=(P'+string(n-1)+'(x))+(z(n+1).*(f'+string(n)+'(x)))');


    w = z(1).*x.^0;
    for i=1:n
    deff('y=flix(x)','y=(f'+string(i)+'(x))');
    w = w + z(i+1).*flix(x)
    end


p = lagr(x,w);
disp(x)
disp(w)
endfunction

function z=graf_pol(a,b,pol,cor)
// Traça curva do polinomio descrito no vetor pol
// de forma que pol = [ a0 a1 ...]  
// representa o polinomio 
// P(x) = a0 + a1*x + ...  
//
// O grafico é gerado baseado numa malha regular de mil
// pontos no intervalo [a b]
n=max(size(pol))-1;
yp=[];   // dados para curva de P
x_=a : (b-a)/1000 : b;
// matriz com potencias de x_
aux=[];
for i = 0:n,
aux = [ aux x_'^i]
end

yp = (pol*aux')'
plot(x_,yp,cor)
z=[]
endfunction




function z=lagr(x,y)
// z=lagr(x,y)
// 
// Cria polinomio de lagrange
// Entradas:  x = pontos a interpolar x = [x0 x1 ... xn] 
//            y = [f(x0) f(x1) ... f(xn)]
// Saida:     
//            z = [ a0 a1 ...] tais que 
//            P(x) = a0 + a1*x + ...  eh o pol interpolador de Lagange
//
// Eduardo. 14/03/10
c=size(x,2);
M=[];
for k=1:c,

 l = 1;
 for j=1:c,
 if j~=k,
   l = inv(x(k)-x(j))*convol([-x(j) 1],l);
 end
 end
M = [ M; y(k)*l];

end

z = sum(M,1);
endfunction
