function [table, tabLogic] = NonComGaussElim(myperms, table, tabLogic)
%
% NONCOMGAUSSELIM: A function used for computing the non-commutive gaussian
% elimination algorithm on a finitely generated subset of the symmetric
% group.
%
% [table,tabLogic] = NonComGaussElim(myperms,table,tabLogic)
%
%  INPUT:  myperms - [cell] A cell whose entries consist of the generators
%                    of the group. The format of these entries should be
%                    given via the "image" notation. Namely, if sigma is
%                    a generator in Sn, then express sigma as
%                    [sigma(1), sigma(2), ..., sigma(n)]
%            table - [3d array, non-user input] The result of the
%                    calculations. Necessary for recursion and needn't be
%                    supplied by the end-user.
%         tabLogic - [2d array, non-user input] The logical array
%                    consisting of the non-identity indices for which
%                    elements have been inserted into the table. Need not
%                    be supplied by the end-user.
% OUTPUT:    table - [3d array] The final result of the calculations.
%                    table(i,j,:) consists of the permutation element
%                    corresponding to pivot i with image j. 
%         tabLogic - [2d array] The logical array
%                    consisting of the non-identity indices for which
%                    elements have been inserted into the table.
%
% Example: Consider the symmetric group S4 and the generators
%          g1 = [2 3 1 4] and g2 = [2 1 4 3]. The code is executed as
%          >>myperms = {[2 3 1 4], [2 1 4 3]};
%          >>[table, tabLogical] = NonComGaussElim(myperms);

% Version : 1.2
% Author  : Tyler Holden
% Date    : September 29, 2011

%    To Do: 
% 1) Add step to pre-process the input permutations. In particular, if
%    the identity element is given as one of the generators, this algorithm
%    will terminate prematurely. This is because the recursive halting
%    criterion is the identity permutation
% 2) Optimize for parallel computing. There are not many places wherein
%    such improvements could be made, since the algorithm is recursive and
%    heavily dependent on neighbouring iterations. However, some
%    improvements may be made to computing products and inverses of
%    elements.

try
    
    % Return condition for recursion.
    if isempty(myperms)
        return
    end %End if

    sizeGrp = length(cell2mat(myperms(1)));

    % Since the user is not expected to input the table, the following
    % condition checks that we are at the top level of the recursion and
    % creates the table as necessary.
    if nargin < 2
        table = zeros(sizeGrp,sizeGrp,sizeGrp);
        tabLogic = zeros(sizeGrp);
        %     for i=1:sizeGrp
        %         table(i,i,:) = 1:sizeGrp;
        %     end %End for(i)
    end %End if

    % Here we pull out the first generator and calculate its pivot point.
    curGen = cell2mat(myperms(1));
    [piv,img] = findPivot(curGen);

    % If curGen is the identity, piv will be empty. In this case, we simply
    % return to the call function. This is used in sublevels of the recursion.
    % However, this will result in improper calculations if the identity is
    % included in the original set of permutations. In a future version of this
    % code, one could "pre-process" the input permutations to remove the
    % identity elements. This would remove the erroneous calculations and the
    % recursion check would stay consistent.
    if isempty(piv)
        return;
    end %End if

    % Pull the permutation from the table. Afterwards, check if it is empty. If
    % so, input curGen into the table. Otherwise, pre-multiply and apply
    % recursively to the singleton permutation.
    sigij(1:sizeGrp) = table(piv,img,1:sizeGrp);

    if isempty(find(sigij))
        table(piv,img,:) = curGen;
        tabLogic(piv,img) = 1;
        % Only when something new has been added to the table do we do the
        % multiplication step. Note that it is important to do multiplication
        % on both sides, otherwise this will not work.
        for it1=sizeGrp:-1:1;
            for it2=sizeGrp:-1:it1
                if tabLogic(it1,it2)
                    sigij(1:sizeGrp) = table(it1,it2,1:sizeGrp);
                    %First multiplication.
                    newPerm = permMult(sigij, curGen);
                    [table,tabLogic] = NonComGaussElim({newPerm}, table, tabLogic);
                    %Second multiplication.
                    newPerm = permMult(curGen, sigij);
                    [table,tabLogic] = NonComGaussElim({newPerm}, table, tabLogic);
                end %End if
            end %End for(it2)
        end %End for (it1)
    else
        newPerm = permMult(permInv(sigij),curGen);
        [table, tabLogic] = NonComGaussElim({newPerm}, table, tabLogic);
    end %end if

    %Once the generator and its multiples have been added, we remove it from
    %the list of permutations and apply the algorithm recursively.
    [table, tabLogic] = NonComGaussElim(myperms(2:end), table, tabLogic);
catch ME
    disp('Error');
end

%-------------------------------Subroutines-------------------------------%

function [pivot, image] = findPivot(permut)
%
% Calculates the pivot of the permutation. If the permutation is identity,
% pivot and image return empty.
pivot = find(permut ~= 1:length(permut),1);
image = permut(pivot);

function prod = permMult(perm1, perm2)
%
% Multiplies two permutations. This is simply the composition of the
% permtuations acting on the set {1,...,n}. That is, if sigma and tau are
% permutations, then sigma*tau acting on the set {1,...,n} is given in
% "image notation" as
% sigma*tau = [sigma(tau(1)), sigma(tau(2)), ... , sigma(tau(n)) ]
prod = perm1(perm2);

function inverse = permInv(perm)
%
% Calculates the inverse element. If sigma is the permutation to be
% inverted, it is written as [sigma(1), sigma(2), ... , sigma(n)]. The
% inverse of this function is just done by placing j in the sigma(j)
% position. 
inverse(perm) = 1:length(perm);