function [p] = MatchLoops(A)
tic
[m,n] = size(A);
Observed = trace(A);
Index = 0; permutation = zeros(1, n);
for i = 1:n-1
    permutation(i) = i;
end
Position = n - 1;
 
while (Position ~= 1) || (permutation(Position) <= n)        
    Position = Position + 1;
    fathom = 0;
    while fathom == 0
        Redundancy = 1;
        while ((Redundancy == 1) && (permutation(Position) <= n))
            permutation(Position) = permutation(Position) + 1;
            Redundancy = 0;
            for k = 1:(Position - 1)
                if permutation(Position) == permutation(k)
                    Redundancy = 1;
                end
            end
        end
        if (Position == 1) && (permutation(Position) > n)
            break;
        end           
        if (Position > 1) && (permutation(Position) > n)
            permutation(Position) = 0;
            Position = Position - 1;
        else
            if Position == n
                stat = 0;
                for i = 1:n
                    stat = stat + A(i, permutation(i));
                end
                if stat >= Observed
                    Index = Index + 1;
                end;
            else
                PREstat = 0; us = 0; greatertest = 0; lesstest = 0;
                for i = 1:n
                    found = false;
                    for j = 1:Position
                        if permutation(j) == i
                            found = true;
                            PREstat = PREstat + A(j, permutation(j));
                        end
                    end
                    if ~found
                        largest = A(Position + 1, i); smallest = A(Position + 1, i);
                        for k = Position + 1:n
                            if largest < A(k, i)
                                largest = A(k, i);
                            end
                            if smallest > A(k, i)
                                smallest = A(k, i);
                            end                                
                        end
                        greatertest = greatertest + smallest;
                        lesstest = lesstest + largest;
                    end                    
                end
                fathom = 1;
                % fathom = GREATERTEST
                if PREstat + greatertest >= Observed
                   fathom = 0; Index = Index + factorial(n - Position);
                else
                    % fathom = LESSTEST
                    if PREstat + lesstest < Observed
                        fathom = 0;
                    end;
                end
            end
        end     % No Retraction
    end     % fathom loop
end     % Termination loop
p = (Index/factorial(n));
Index=Index
toc