function [X, flag, K] = edm(D)
%       [X, flag] = edm(D)
%   
%   INPUT:
%   D: matriz de distancias (Dij = dij^2) completa
%
%   OUTPUT:
%   flag: 0 se D nao e' EDM e 1 caso contrario
%   X: se D e' EDM, X retorna uma realizacao na dimensao de realizacao K
%
    eps1 = 1e-8;
    n = length(D);
    
    % matriz de centragem
    J = eye(n) - (1/n)*ones(n,n);
    
    % transformacao linear
    G = -(1/2)*J*D*J;
    
    % decomposicao espectral de G
    [Q,La] = eig(G,'vector');
    
    if min(La) < -eps1
        % D nao e' EDM
        flag=0;
        X = [];
        return
    end
    
    [Lo, idx] = sort(La,'descend');
    K = length(find(Lo > eps1));
    X = diag(sqrt(Lo(1:K)))*Q(:,idx(1:K))';

    flag = 1; % D e' EDM

end