McGill University
School of Computer Science

Computer Science COMP 199 (Winter term)
Excursions in Computing Science

%MATLABpak08cI: groupgenMatrix.m and supply progs, groupgen.m and supply progs.
%
%	groupgenMatrix.m		insertMatColl.m
%	insertMatColl.m			equalmat.m
%
%	groupgen.m			prodcyc.m
%	groupgen.m			invcyc.m
%	groupgen.m			insertgroup.m
%	prodcyc.m			array2cyc.m
%	array2cyc.m			cyc2array.m
%	cyc2array.m			commacycles.m
%	array2cyc.m			prodperm.m
%	array2cyc.m			invperm.m
%	invperm.m			array2cyc.m
%	invperm.m			cyc2array.m
%	invperm.m			invperm.m
%	prodcyc.m			cyc2array.m
%	prodcyc.m			prodperm.m
%	insertgroup.m			equalcyc.m

% function group = groupgenMatrix(generator)		THM		071102
% Finds closure of identity element and group elements in  generator
% E.g., generator = {[0,-1;1,0]}; rot4 = groupgenMatrix(generator); gives cell
%  array rot4: {[1,0;0,1],[0,-1;1,0],[0,1;-1,0],[0,-1;-1,0]}
% Uses brute force and
%   matColl = insertMatColl(element,matColl)
function group = groupgenMatrix(generator)
ng = size(generator); ng = ng(1);
sizmat = size(generator{1});	% all elements are square matrices of same size
group{1} = diag(diag(ones(sizmat(1)),0),0);
			% unit matrix same size as generator mat
% insert the generator elements and their inverses
for k = 1:ng
  group = insertMatColl(generator{k},group);
  group = insertMatColl(generator{k}^-1,group);
end %for k = 1:ng
% repeat very inefficiently until group stops growing
sg1 = size(group);
while true
  sg = sg1;
  for k = 1:sg(2)
    for m = 1:sg(2)
      group = insertMatColl(group{k}*group{m},group);
      group = insertMatColl(group{m}*group{k},group);
    end %for m = 1:sg(2)
  end %for k = 1:sg(2)
  sg1 = size(group)		% this takes a long time, so report size so far
%	for k = 1: sg1(2)
%	group{k}						% test only
%	end %for k = 1: sg1{2}
  if sg1(2) == sg(2), break, end
end %while true

% function MatColl = insertMatColl(element,MatColl)		THM	071102
% Uses
%   eqmat = equalmat(opmatrix1,opmatrix2)
% to find if the  element  is already in the cell  MatColl  and inserts if not
function MatColl = insertMatColl(element,MatColl)
sg = size(MatColl);
there = false;
for k = 1:sg(2)
  there = there | equalmat(element,MatColl{k});
  if there, break, end
end %for k = 1:sg(2)
if ~there
  MatColl{sg(2)+1} = element;	% insert
end %if ~there

% function eqmat = equalmat(opmatrix1,opmatrix2)	THM		071102
% eqmat = opmatrix1 == opmatrix2
function eqmat = equalmat(opmatrix1,opmatrix2)
% not equal if sizes differ
n1 = size(opmatrix1);
n2 = size(opmatrix2);
eqmat = n1(1)==n2(1) & n1(2)==n2(2);
if eqmat
% compare
% compare = opmatrix1 == opmatrix2;	% produces matrix of 1s, 0s
  for k = 1:n1(1)
    for m = 1:n1(2)
      if abs(opmatrix1(k,m))<10^-6
        eqmat = eqmat & abs(opmatrix1(k,m)-opmatrix2(k,m))<10^-6; else
      eqmat = eqmat & abs((opmatrix1(k,m)-opmatrix2(k,m))/opmatrix1(k,m))<10^-6;
      end %if opmatrix1(k,m)<10^-6
      if ~eqmat, break, end
    end %for m = 1:n1(2)
  end %for k = 1:n1(1)
end %if eqmat


% function group = groupgen(n,generator)		THM		071029
% Finds closure of identity element and group elements in  generator
% E.g., generator = {'(12345)','(123)','(12)(34)'}  NB or '(1,2,3,4,5)', ..
% Uses brute force and
%   opcycles = prodcyc(n,opcycles2,opcycles1)
%   opcycles = invcyc(n,opcycles1)
%   group = insertgroup(n,element,group1)
% n is the highest integer in the cycle notation: permutations of integers 1--n
function group = groupgen(n,generator)
ng = size(generator); ng = ng(2);
group{1} = '()';
% insert the generator elements and their inverses
for k = 1:ng
  group = insertgroup(n,generator{k},group);
  group = insertgroup(n,invcyc(n,generator{k}),group);
end %for k = 1:ng
% repeat very inefficiently until group stops growing
sg1 = size(group);
while true
  sg = sg1;
  for k = 1:sg(2)
    for m = 1:sg(2)
      group = insertgroup(n,prodcyc(n,group{k},group{m}),group);
      group = insertgroup(n,prodcyc(n,group{m},group{k}),group);
    end %for m = 1:sg(2)
  end %for k = 1:sg(2)
  sg1 = size(group)		% this takes a long time, so report size so far
  if sg1(2) == sg(2), break, end
end %while true

% function opcycles = prodcyc(n,opcycles2,opcycles1)	THM		071029
% Uses
%   opcycles = array2cyc(oparray)
%   oparray = cyc2array(n,opcycles)
%   oparray = prodperm(oparray2,oparray1)
% to find products of permutations expressed as cycles
function opcycles = prodcyc(n,opcycles2,opcycles1)
opcycles = array2cyc(prodperm(cyc2array(n,opcycles2),cyc2array(n,opcycles1)));

% function opcycles = array2cyc(oparray)	THM			071026
% converts array notation for permutations to cycles
% e.g., 6:(2,3,4,5) <-	1 2 3 4 5 6 
%			1 3 4 5 2 6	:oparray(1,1:6) the permutation
%			0 1 1 1 1 0	:oparray(2,1:6) the cycles
% e.g., 6:(2,5)(3,4) <-	1 2 3 4 5 6 
%			1 5 4 3 2 6	:oparray(1,1:6) the permutation
%			0 1 2 2 1 0	:oparray(2,1:6) the cycles
% Note that the cycles row is useful as a flag workspace. It is kept anyway
% and internally regenerated in canonical order of a) least term starts cycle
% and b) cycles ordered by their least term.
% companion functions:
%   oparray = cyc2array(n,opcycles)
%   oparray = prodperm(oparray2,oparray1)
%   oparray = invperm(oparray1)
% do a pass for each cycle
function opcycles = array2cyc(oparray)
soa = size(oparray); n = soa(2);
% initialize cycle labels to 0
for k = 1:n; oparray(2,k) = 0; end
% output the cycles and rehuild them in oparray in canonical order
% (note that these get lost since oparray is not output)
cyc = 1;
opcycles = '';
for k = 1:n
  if oparray(2,k)==0		% use cycle label to flag if visited
    k1 = oparray(1,k);
    if k1 ~= k			% ignore singleton cycles
      opcycles = [opcycles '(' int2str(k)];
      oparray(2,k) = cyc;
      while true
	opcycles = [opcycles ',' int2str(k1)];
	oparray(2,k1) = cyc;
	k1 = oparray(1,k1);	% next in cycle
	if k1==k, break, end
      end %while true
      opcycles = [opcycles ')'];
      cyc = cyc + 1;
    end %if k1 ~= k
  end %if oparray(2,k)==0
end %for k = 1:n
if isempty(opcycles), opcycles = '()'; end	% identity

% function oparray = cyc2array(n,opcycles)	THM			071024
% converts cycle notation for permutations to arrays
% e.g., 6,(2,3,4,5) ->	1 2 3 4 5 6 
%			1 3 4 5 2 6	:oparray(1,1:6) the permutation
%			0 1 1 1 1 0	:oparray(2,1:6) the cycles
% e.g., 6,(2,5)(3,4) ->	1 2 3 4 5 6 
%			1 5 4 3 2 6	:oparray(1,1:6) the permutation
%			0 1 2 2 1 0	:oparray(2,1:6) the cycles
% Note that the cycles row is useful as a flag workspace. It is kept anyway
% and internally regenerated in canonical order of a) least term starts cycle
% and b) cycles ordered by their least term.  cyc2array()  records cycles as
% read from  opcycles
% uses
%   opcycles = commacycles(n, opcycles);
% companion functions:
%   opcycles = array2cyc(oparray)
%   oparray = prodperm(oparray2,oparray1)
%   oparray = invperm(oparray1)
% pseudocode based on a 7-state automaton (states t,u,v,w,x,y,z)
% (this algorithm ignores noise in input by staying in current state except
% for states v, y, which escape to states w, z, respectively)
%	state <- t	initial state
%	cyc <- 0	cycle identifier, to be incremented from 0
%	oparray <- [1,2,3,4,5,6,..;0,0,0,0,0,0,..]	initialized
%	[current state, input char, next state]
%	t ( u	cyc++	num <- 0    num is the numeric between ( or , and , or )
%	u ) t	cyc--			identity is ()
%	u d v	num <- num*10 + d	d is the digit read
%	v d v	num <- num*10 + d	d is the digit read
%	w d v	num <- num*10 + d	d is the digit read
%				(get to states w, z only if noise in input)
%	v ) t	cyc--	do nothing else: singleton cycles already initialized as
%			cycle 0   (if no initialization: oparray(1,num) <- num
%							 oparray(2,num) <- 0 )
%	v , x	indx <- num	this is the location the next num will replace
%		num1st <- num	this number will replace the final loc in cycle
%		num <- 0	initialize for next number
%	y d y	num <- num*10 + d	d is the digit read
%	z d y	num <- num*10 + d	d is the digit read
%	y d x	num <- num*10 + d	d is the digit read
%				(get to states w, z only if noise in input)
%	y , x	oparray(1,indx) <- num
%		oparray(2,indx) <- cyc
%		indx <- num	num <- 0
%	y , t	oparray(1,indx) <- num
%		oparray(2,indx) <- cyc
%		indx <- num	num <- 0
%		oparray(1,indx) <- num1st
%		oparray(2,indx) <- cyc
function oparray = cyc2array(n,opcycles)
% expand commaless cycles on 10 or fewer items
opcycles = commacycles(n, opcycles);
% initialize
t = 0; u = 1; v = 2; w = 3; x = 4; y = 5; z = 6;
state = t; cyc = 0;
for k = 1:n; oparray(1,k) = k; oparray(2,k) = 0; end
% run automaton
soc = size(opcycles);
for k = 1:soc(2)
  nextsymb = opcycles(k);			%; after test
  switch state
  case t
    switch nextsymb
    case '('
      cyc = cyc + 1;
      num = 0;
      state = u;				%; after test
   %otherwise noise: do nothing
    end %switch nextsymb
  case u
    if nextsymb==')'
      state = t;				%; after test
      cyc = cyc - 1;				%; after test
    elseif '0' <= nextsymb & nextsymb <= '9'
      digit = str2num(nextsymb);
      num = num*10 + digit;			%; after test
      state = v;				%; after test
   %else noise: do nothing
    end %if nextsymb==')' elseif '0' <= nextsymb & nextsymb <= '9'
  case v
    switch nextsymb
    case ')'
      cyc = cyc - 1;
      % singleton cycle, taken care of by initialization
      state = t;				%; after test
    case ','
      indx = num;	%this is the location the next num will replace
      num1st = num;	%this number will replace the final loc in cycle
      num = 0;		%initialize for next number
      state = x;				%; after test
    otherwise
      if '0' <= nextsymb & nextsymb <= '9'
	digit = str2num(nextsymb);
	num = num*10 + digit;			%; after test
	state = v;				%; after test
      else %noise:
        state = w;				%; after test
      end %if '0' <= nextsymb & nextsymb <= '9'
    end %switch nextsymb
  case w %noise state
    if '0' <= nextsymb & nextsymb <= '9'
      digit = str2num(nextsymb);
      num = num*10 + digit;			%; after test
      state = v;				%; after test
   %else more noise: do nothing
    end %if '0' <= nextsymb & nextsymb <= '9'
  case x
    if '0' <= nextsymb & nextsymb <= '9'
      digit = str2num(nextsymb);
      num = num*10 + digit;			%; after test
      state = y;				%; after test
   %else noise: do nothing
    end %if '0' <= nextsymb & nextsymb <= '9'
  case y
    switch nextsymb
    case ')'
      oparray(1,indx) = num;
      oparray(2,indx) = cyc;
      indx = num; num = 0;
      oparray(1,indx) = num1st;
      oparray(2,indx) = cyc;
      state = t;				%; after test
    case ','
      oparray(1,indx) = num;
      oparray(2,indx) = cyc;
      indx = num;	%this is the location the next num will replace
      num = 0;		%initialize for next number
      state = x;				%; after test
    otherwise
      if '0' <= nextsymb & nextsymb <= '9'
	digit = str2num(nextsymb);
	num = num*10 + digit;			%; after test
	state = y;				%; after test
      else %noise:
        state = z;				%; after test
      end %if '0' <= nextsymb & nextsymb <= '9'
    end %switch nextsymb
  case z
    if '0' <= nextsymb & nextsymb <= '9'
      digit = str2num(nextsymb);
      num = num*10 + digit;			%; after test
      state = y;					%; after test
   %else more noise: do nothing
    end %if '0' <= nextsymb & nextsymb <= '9'
  end %switch state
end %for k = 1:soc(2)  i.e. size(opcycles)(2)

% function opcycles = commacycles(n, opcycles1)		THM		071029
% Expands opcycles with no commas for n<=10 by adding them
% pseudocode based on a 2-state automaton (states y,z)
% (ignores all but ',' and digits by not changing state; copies all input)
% if n <= 0
%   state = y
%   [current state, input char, next state]
%   y d z
%   z d z	insert ','
%   z , y
% Used by
%  oparray = cyc2array(n,opcycles)
%  sims = allsims(n,element,group)
function opcycles = commacycles(n, opcycles1)
if n > 10
  opcycles = opcycles1
else
  soc = size(opcycles1);
  y = 0; z = 1;
  state = y;
  k = 1;
  for k1 = 1:soc(2)
    nextsymb = opcycles1(k1);
    switch state
      case y
	if '0'<= nextsymb & nextsymb <= '9'
	  state = z;
	end %if '0'<= nextsymb & nextsymb <= '9'
      case z
	if nextsymb == ','
	  state = y;
	elseif '0'<= nextsymb & nextsymb <= '9'
	  opcycles(k) = ',';
	  k = k + 1;
	  state = z;
	end %if '0'<= nextsymb & nextsymb <= '9'
    end %switch state
    opcycles(k) = opcycles1(k1);	% copy input
    k = k + 1;
  end %for k1 = 1:soc(2)
end %if n > 10

% function oparray = prodperm(oparray2,oparray1)	THM		071026
% combines two permutations ordered oparray1 then oparray2
% E.g. 123456 then 123456 gives 123456 and we recreate the cycles 123456
%      154326      134526       125436                            125436
%                                                                 001010
% Note that the cycles row is useful as a flag workspace. It is kept anyway
% and internally regenerated in canonical order of a) least term starts cycle
% and b) cycles ordered by their least term.
% companion functions:
%   oparray = cyc2array(n,opcycles)
%   opcycles = array2cyc(oparray)
%   oparray = invperm(oparray1)
function oparray = prodperm(oparray2,oparray1)
soa = size(oparray1); n = soa(2);	% assume soa == size(oparray2)
% product
for k = 1:n
  oparray(1,k) = oparray2(1,oparray1(1,k));
end %for k = 1:n
% initialize cycle labels to 0
for k = 1:n; oparray(2,k) = 0; end
% rebuild cycles in canonical order
cyc = 1;
for k = 1:n
  if oparray(2,k)==0		% use cycle label to flag if visited
    k1 = oparray(1,k);
    if k1 ~= k			% ignore singleton cycles
      oparray(2,k) = cyc;
      while true
	oparray(2,k1) = cyc;
	k1 = oparray(1,k1);	% next in cycle
	if k1==k, break, end
      end %while true
      cyc = cyc + 1;
    end %if k1 ~= k
  end %if oparray(2,k)==0
end %for k = 1:n

% function oparray = invperm(oparray1)		THM			071026
% E.g. 123456 inverts to 123456: just swap rows 1, 2 and resort (but it's
%      156243            146523        easier to do than that)
%      021221            021221
% Note that the cycles row is useful as a flag workspace. It is kept anyway
% and internally regenerated in canonical order of a) least term starts cycle
% and b) cycles ordered by their least term.  invperm() leaves it alone
% companion functions:
%   oparray = cyc2array(n,opcycles)
%   opcycles = array2cyc(oparray)
%   oparray = prodperm(oparray2,oparray1)
function oparray = invperm(oparray1)
soa = size(oparray1); n = soa(2);
for k = 1: n
  oparray(1,oparray1(1,k)) = k;
  oparray(2,oparray1(1,k)) = oparray1(2,k);
end %for k = 1: n

% function opcycles = invperm(n,opcycles1)	THM			071029
% Uses
%   opcycles = array2cyc(oparray)
%   oparray = cyc2array(n,opcycles)
%   oparray = invperm(oparray1)
% to find inverse of a permutation expressed as cycles
function opcycles = invperm(n,opcycles1)
opcycles = array2cyc(invperm(cyc2array(n,opcycles1)));

% function group = insertgroup(n,element,group)		THM		071029
% Uses
%   eqcyc = equalcyc(n,opcycles1,opcycles2)
% to find if the  element  is already in the cell  group  and insert if not
function group = insertgroup(n,element,group)
sg = size(group);
there = false;
for k = 1:sg(2)
  there = there | equalcyc(n,element,group{k});
  if there, break, end
end %for k = 1:sg(2)
if ~there
  group{sg(2)+1} = element;
end %if ~there
`
% function eqcyc = equalcyc(n,opcycles1,opcycles2)	THM		071029
% eqcyc = opcycles1 == opcycles2
% NB first put both into canonical form using
%   oparray = cyc2array(n,opcycles)
%   opcycles = array2cyc(oparray)
function eqcyc = equalcyc(n,opcycles1,opcycles2)
% put into canonical form
opcycles1 = array2cyc(cyc2array(n,opcycles1));
opcycles2 = array2cyc(cyc2array(n,opcycles2));
% not equal if sizes differ
n1 = size(opcycles1); n1 = n1(2);
n2 = size(opcycles2); n2 = n2(2);
eqcyc = n1==n2;
if eqcyc
% compare
  compare = opcycles1 == opcycles2;	% produces array of 1s, 0s
  for k = 1:n1
    eqcyc = eqcyc & compare(k);		% must all be 1s
    if ~eqcyc, break, end
  end %for k = 1:n1
end %if eqcyc