function [L, U, p, sig] = ludecompsig(A)
%LUDECOMPSIG LU factorization
% [L, U, p, sig] = ludecompsig(A) produces a unit lower triangular
% matrix L, an upper triangular matrix U, and a
% permutation vector p and a scalar sig, so that L*U = A(p,:) and
% sig = +1 or -1 if p is an even or odd permutation.

    [~, n] = size(A);
    p = (1:n)';
	sig = 1;

    for k = 1:n-1
    
        % Find the index of the largest element below diagonal in k-th column
        [~, m] = max(abs(A(k:n, k)));
        m = m + k - 1;
    
        % Swap pivot row
        if (m ~= k)
            A([k m], :) = A([m k], :);
            p([k m]) = p([m k]);
            sig = -sig;  % a swap changes the parity of permutation
        end
    
        i = k+1:n;
		% Compute multipliers
        A(i, k) = A(i, k)/A(k,k);
		% Update the remainder of the matrix:
		%   the operation A(i,k)*A(k,j) multiplies a column vector 
		%   by a row vector to produce a square matrix of 
		%   order n-k; this matrix is then subtracted from the
		%   submatrix of the same size in the bottom right corner of A
        j = k+1:n;
        A(i, j) = A(i, j) - A(i,k)*A(k,j);

    end


    % Separate results
    U = triu(A);
    L = tril(A);
    L(1:n+1:end) = 1.;    

end
