function [phim,phiv,EC,alphac,GC]=cholesky_full(xc,N,p)
%  cholesky decomposition
%
%  first compute phim(1,1)...phim(p,p)
%  next compute phiv(1,0)...phiv(p,0)
%
    for i=1:p
        for k=1:p
            phim(i,k)=sum(xc(p+1-i:p+N-i).*xc(p+1-k:p+N-k));
        end
    end
    
    for i=1:p
        phiv(i)=sum(xc(p+1-i:p+N-i).*xc(p+1:p+N));
        phiz(i)=sum(xc(p+1:p+N).*xc(p+1-i:p+N-i));
    end
    phi0=sum(xc(p+1:p+N).^2);
%
% cholesky decomposition
% solve matrix equation phim=v d v' where v is lower triangular and d is
% diagonal matrix (represented as a row vector)
    v=diag(ones(1,p),0);
%
% special case for j=1
    d(1)=phim(1,1);
    v(2:p,1)=phim(2:p,1)/d(1);
    for j=2:p-1
        i=j;
        vsum=sum(v(i,1:i-1).^2.*d(1:i-1));
        d(i)=phim(i,i)-vsum;
        for i=j+1:p
            vsum=sum(v(i,1:j-1).*d(1:j-1).*v(j,1:j-1));
            v(i,j)=(phim(i,j)-vsum)/d(j);
        end
    end
    j=p;
    i=j;
    vsum=sum(v(i,1:j-1).^2.*d(1:j-1));
    d(i)=phim(i,i)-vsum;
    %v
    %d
%
% unwrap solution from v and d
% solve matrix equation v d v' alpha = psi
    y(1)=phiv(1);
    for i=2:p
        vsum=sum(v(i,1:i-1).*y(1:i-1));
        y(i)=phiv(i)-vsum;
    end
    i=p;
    alphac(i)=y(i)/d(i);
    for i=p-1:-1:1
        vsum=sum(v(i+1:p,i).*alphac(i+1:p)');
        alphac(i)=y(i)/d(i)-vsum;
    end
    alphac=alphac';           
    phiv=phiv';
    EC=phi0-sum(alphac'.*phiz);
    GC=sqrt(EC);