function Y = matrix_sign(X)
% Y = MATRIX_SIGN(X) computes the matrix sign function by means of the
% scaled Newton's method
% X: a square matrix
% Y: the sign of X
n = size (X,1);
tol = 1e-13;
kmax = 20;
err = 1;
Y = X;
k = 0;
while err>tol && k<kmax
    [L,U,P] = lu(Y);
    c = abs(prod(diag(U)))^(-1/n);
    Y = Y*c;
    Z = L\P;
    Z = (c*U)\Z;
    err = norm(Y-Z,1)/norm(Y,1);
    Y = (Y+Z)*0.5;
    k = k + 1;
end
if k == kmax
    disp('Warning: reached the maximum number of iterations')
end