clear;

function r=pscalL2(F,G)
    E=F*G'
    N=length(F);
    r=0;
    for i=1:N
        for j=1:N
            k=(i-1)+(j-1)
            r=r+E(i,j)*(1-(-1)^(k+1))/(k+1);
        end
    end
endfunction

function P=orthoL2can(N)
    P=diag(ones(N,1));
    P(:,1)=P(:,1)/sqrt(pscalL2(P(:,1),P(:,1)));
    for i=2:N
        for j=1:i-1
            P(:,i)=P(:,i)-pscalL2(P(:,i),P(:,j))*P(:,j);
        end
        P(:,i)=P(:,i)/sqrt(pscalL2(P(:,i),P(:,i)))
    end
endfunction



function P=Legendre(N)
    P=diag(ones(N,1));
    for i=2:N-1
        P(:,i+1)=((2*(i-1)+1)*[0;P(1:N-1,i)] -(i-1)*P(:,i-1))/i;
    end 
    for i=1:N
        P(:,i)=P(:,i)/sqrt(pscalL2(P(:,i),P(:,i)))
    end

endfunction

function Y=evaluatepoly(P,X)
    N=length(P)
    Y=zeros(1,length(X));
    for i=1:N
        Y=Y+P(i)*X^(i-1);
    end

endfunction

function A=approxL2(N)
    A=zeros(N,1);
    P=Legendre(N);
    for i=1:N
        pscal=0; 
        for j=1:N
            if(modulo(j-1,2)==0)
                pscal=pscal+(4/(4*(j-1)+1))*P(j,i);
            end
        end
        A=A+pscal*P(:,i);
    end

endfunction



disp(pscalL2([1 0 0]',[0 0 1]'));
//on trouve 0.6666667
disp(pscalL2([1 1 0 0]',[0 1 0 1]'));//(1+x), (x+x^3)
// on trouve 1.0666667  , i.e. 2/3+2/5
N=5

PP=Legendre(N)
P=orthoL2can(N)

erreur=norm(P-PP,'inf')
X=linspace(-1,1,100);

clf
for i=1:N
    //plot(X,evaluatepoly(PP(:,i),X));
    //plot(X,evaluatepoly(P(:,i),X));
    plot2d(X,evaluatepoly(approxL2(i),X),style=1);
end

plot(X,sqrt(abs(X)).^(-1),'-')

