#This code creates the hyperoctahedral group or signed permutation group, note that if you want B_n then you will get a group on 2n letters (replace n+1, n+2 with \bar{1} \bar{2} etc #
#Example 
#gap>G:= CubeGroup(3);
#Group([ (1,4), (1,2)(4,5), (2,3)(5,6) ])
#To get the elements of this group use 
#Elements(G);
#[ (), (3,6), (2,3)(5,6), (2,3,5,6), (2,5), (2,5)(3,6), (2,6,5,3), (2,6)(3,5), (1,2)(4,5), (1,2)(3,6)(4,5), 
 # (1,2,3)(4,5,6), (1,2,3,4,5,6), (1,2,4,5), (1,2,4,5)(3,6), (1,2,6,4,5,3), (1,2,6)(3,4,5), (1,3,2)(4,6,5), (1,3,5,4,6,2), 
 # (1,3)(4,6), (1,3,4,6), (1,3,2,4,6,5), (1,3,5)(2,4,6), (1,3)(2,5)(4,6), (1,3,4,6)(2,5), (1,4), (1,4)(3,6), 
 # (1,4)(2,3)(5,6), (1,4)(2,3,5,6), (1,4)(2,5), (1,4)(2,5)(3,6), (1,4)(2,6,5,3), (1,4)(2,6)(3,5), (1,5,4,2), 
 # (1,5,4,2)(3,6), (1,5,6,4,2,3), (1,5,6)(2,3,4), (1,5)(2,4), (1,5)(2,4)(3,6), (1,5,3)(2,6,4), (1,5,3,4,2,6), 
 # (1,6,5,4,3,2), (1,6,2)(3,5,4), (1,6,4,3), (1,6)(3,4), (1,6,5)(2,4,3), (1,6,2,4,3,5), (1,6,4,3)(2,5), (1,6)(2,5)(3,4) ]
#gap> 
#Which returns all 48 elements of B_3
CubeGroup:= function(n)
   local s, i;
  s:= List([1..n], x->[1..2*n]);
     s[1]{[1, n+1]}:= [n+1, 1];
   for i in [2..n] do
s[i]{[i-1, i, n+i-1, n+i]}:= [i, i-1, n+i, n+i-1]; 
     od;
     return Group(List(s, PermList), ());
   end;

#Just like your classic conj classes but only conjugates with elements from H#
# Example use as follows, if G:=SymmetricGroup(3) and H:=Stabilizer(G,[3],OnSets), then H is the Young Subgroup of G isomorphic to S_2 x S_1
#gap> HConjClass(G,H);
#[ [ () ], [ (2,3), (1,3) ], [ (1,2) ], [ (1,2,3), (1,3,2) ] ]

HConjClass:=function(G,H)
local i, j, k, a, b, c, d, p,Y, y, x, X;
i:=1;
j:=1;
b:=[];
c:=[];
y:= Size(G);
Y:= Elements(G);
x:=Size(H);
X:= Elements(H);

while i <=y do
	j:=1;
	#Print("i= ",i,"\n");
	b:=[];
	k:=1;
	p:=0;
	while k<=Length(c) do
		if Y[i] in c[k] then
			p:=1;
			k:=Length(c)+1;
		else
			p:=0;
		fi;
		k:=k+1;
	od;
	if p=0 then
		Add(b,Y[i]);
		while j<=x do
			a:=(X[j]^(-1))*Y[i]*(X[j]);
			if a in b then
				Print("");
			else
				Add(b,a);
			fi;
			j:=j+1;
		od;
		if b in c then
			Print("");
		else
			Add(c,b);
		fi;
	else
		Print("");
	fi;
	i:=i+1;
od;

return c;	
end;




#This function is used in Dixons algorithm
#Take a matrix M and Convert all the integers to elements in Z_p#
#gap> M:=[[1,2,3],[0,1,1],[4,5,6]];
#[ [ 1, 2, 3 ], [ 0, 1, 1 ], [ 4, 5, 6 ] ]
#gap> A:=MatToZ(M,3);
#[ [ Z(3)^0, Z(3), 0*Z(3) ], [ 0*Z(3), Z(3)^0, Z(3)^0 ], [ Z(3)^0, Z(3), 0*Z(3) ] ]
#gap> ZToMat(A,3);   
#[ [ 1, 2, 0 ], [ 0, 1, 1 ], [ 1, 2, 0 ] ]



MatToZ:=function(M,p)
local i, j, k;
i:=1;
j:=1;
while j<=Size(M) do
	i:=1;
	while i<=Size(M[j]) do
		M[j][i]:= M[j][i]*Z(p)^0;
		i:=i+1;
	od;
j:=j+1;
od;
return M;
end;

#Take a matrix M and convert all the elements in Z_p to integer#

ZToMat:=function(M,p)
local i, j, k, B;
B:=[];
i:=1;
j:=1;
while j<=Size(M) do
	i:=1;
	B[j]:=[];
	while i<=Size(M[j]) do
		B[j][i]:=IntFFE(M[j][i]);
		i:=i+1;
	od;
j:=j+1;
od;
return B;
end;


#Used in Dixons algorithm
#Here we use gap to find the class multiplication coeffiecients which #
#Gap computes for us by using the character table#
#gap> ClassCoeffsCharTab(SymmetricGroup(2));
#[ [ [ 1, 0 ], [ 0, 1 ] ], [ [ 0, 1 ], [ 1, 0 ] ] ]


ClassCoeffsCharTab:=function(G)
local i, j, k, c, p, h, H, P;
c:=CharacterTable(G);
i:=1;
H:=[];
H:=SortedConjClass(G,c);;
h:=Length(H);
P:=[];
while i<=h do
	j:=1;
	P[i]:=[];
	while j<=h do
		k:=1;
		P[i][j]:=[];
		while k<=h do
			p:=ClassMultiplicationCoefficient(c,i,j,k);
			#Print(" ",p," ");
			P[i][j][k]:=p;
		k:=k+1;
		od;
	j:=j+1;	
	od;
i:=i+1;	
od;
 return P;
end;

#Given a list of conjugacy classes (as a list of Elements) this returns power to n  of these conjugacy classes#
#Example
#gap> a:=[ (3,6), (2,5), (1,4) ]
#gap> gap> PowerConjClass([a],3);
#[ [ (3,6), (2,5), (1,4) ] ];
#gap> PowerConjClass([a],2);
#[ [ () ] ]

PowerConjClass:=function(K,n)
local i, j, I;
i:=1;
I:=[];
j:=1;
while j<=Size(K) do;
	i:=1;
	I[j]:=[];
	while i<=Size(K[j]) do;
		I[j][i]:=K[j][i]^n;
		i:=i+1;
	od;
	I[j]:=Set(I[j]);
j:=j+1;
od;
return I;
end;


#If K_1, K_2, ..., K_k is a list of conjugacy classes then this returns the list [1,2, ..., k'] which corresponds to the list of power conjugacy classes #
#Example 
#G:=CubeGroup(2);
#K:=HConjClass(G,G);
#[ [ () ], [ (2,4), (1,3) ], [ (1,2)(3,4), (1,4)(2,3) ], [ (1,2,3,4), (1,4,3,2) ], [ (1,3)(2,4) ] ]
#gap> PowerConjClassDash(K,2);
#[ 1, 1, 1, 5, 1 ]
# Interpret this as (K_2)^2 = K_1 etc
# gap> InvConjjdash(K);
#gap>[ 1, 2, 3, 4, 5 ]
PowerConjClassDash:=function(K,n)
local i, j, jdash, I, k, r;
jdash:=[];
I:=[];
k:=Size(K);
I:=PowerConjClass(K,n);
i:=1;
j:=1;
while j<=k do;
	i:=1;
	while i<=k do;
		if  IsEqualSet(I[j],K[i]) then
			jdash[j]:=i;
			i:=k;
		fi;
		i:=i+1;
	od;
j:=j+1;
od;
return jdash;
end;


#If K_1, K_2, ..., K_k is a list of conjugacy classes then this returns the list [1,2, ..., k'] which corresponds to the list of inverse conjugacy classes #
#See above for example
InvConjjdash:=function(K)
local i, j, jdash, I, k, r;
jdash:=[];
I:=[];
k:=Size(K);
I:=InverseConjclass(K);
i:=1;
j:=1;
while j<=k do;
	i:=1;
	while i<=k do;
		r:=ClassEqual([I[j]],[K[i]]);
		if r=1 then
			jdash[i]:=j;
			i:=k+1;
		else
			i:=i+1;
		fi;
	od;
j:=j+1;
od;
return jdash;
end;


#THis takes a Group and computes the Class coefficients the long way#
#it should return the same answer as ClassCoeffsCharTab:=function(G)
#gap> G:=CubeGroup(2);
#gap> K:=HConjClass(G,G);
#[ [ () ], [ (2,4), (1,3) ], [ (1,2)(3,4), (1,4)(2,3) ], [ (1,2,3,4), (1,4,3,2) ], [ (1,3)(2,4) ] ]
#gap> ClassCoeffs(K);                
#[ [ [ 1, 0, 0, 0, 0 ], [ 0, 1, 0, 0, 0 ], [ 0, 0, 1, 0, 0 ], [ 0, 0, 0, 1, 0 ], [ 0, 0, 0, 0, 1 ] ], 
#  [ [ 0, 1, 0, 0, 0 ], [ 2, 0, 0, 0, 2 ], [ 0, 0, 0, 2, 0 ], [ 0, 0, 2, 0, 0 ], [ 0, 1, 0, 0, 0 ] ], 
#  [ [ 0, 0, 1, 0, 0 ], [ 0, 0, 0, 2, 0 ], [ 2, 0, 0, 0, 2 ], [ 0, 2, 0, 0, 0 ], [ 0, 0, 1, 0, 0 ] ], 
#  [ [ 0, 0, 0, 1, 0 ], [ 0, 0, 2, 0, 0 ], [ 0, 2, 0, 0, 0 ], [ 2, 0, 0, 0, 2 ], [ 0, 0, 0, 1, 0 ] ], 
#  [ [ 0, 0, 0, 0, 1 ], [ 0, 1, 0, 0, 0 ], [ 0, 0, 1, 0, 0 ], [ 0, 0, 0, 1, 0 ], [ 1, 0, 0, 0, 0 ] ] ]


ClassCoeffs:=function(K)
local i, j, k, x, y, z,l , m, n, h, r, s, t, c, p , q , g, C, M, e, pq, eig, V, B, d, D, Grp, Diag;
h:=[];
c:=[];
M:=[];
k:=Size(K);
#Print("k = ",k,"\n");
#Find the sizes of the H-conjugacy classes#
i:=1;
while i <=k do
	h[i]:=Size(K[i]);
	i:=i+1;
od;

p:=1;
q:=1;
r:=1;
#Here we create set all the class function coefficients to zero#
while p<=k do
	c[p]:=[];
	q:=1;
	while q<=k do
		c[p][q]:=[];
		r:=1;
		while r<=k do
			c[p][q][r]:=[0];
			r:=r+1;
		od;
	q:=q+1;
	od;
p:=p+1;
od;
		
			
#Print(h);#
#We loop through the Conjugacy Classes$
i:=1;
t:=1;
#Print("k =",k,"\n");#
while t <= k do
	j:=1;
	while j<=1 do
		z:=K[t][j];
		r:=1;
		while r<=k do
			l:=1;
			while l<=h[r] do	
				x:=K[r][l];
				s:=1;
				while s<=k do
					m:=1;
					while m<=h[s] do	
						y:=K[s][m];
						if x*y = z then
							c[r][s][t]:=c[r][s][t]+1;
						fi;
					m:=m+1;
					od;
				s:=s+1;
				od;
			l:=l+1;
			od;
		r:=r+1;
		od;		
	j:=j+1;
	od;
t:=t+1;
od;

#Here we pur the c_rst in matrices M_r#
M:=[];
r:=1;
while r<=k do
	M[r]:=[];
	s:=1;
	while s<=k do
		M[r][s]:=[];
		t:=1;
		while t<=k do
			Add(M[r][s],c[r][s][t][1]);
			t:=t+1;
		od;
	s:=s+1;
	od;
r:=r+1;
od;	


return M;
end;

Inverse2:=function(a,p)
local i, j, k;
k:=Size(a);
i:=1;
while i<=Length(a) do
	if a[i]=0*Z(p) then
	Print("");

	else 
	j:=Inverse(a[i]);
	i:=k;
	fi;
i:=i+1;
od;
return j;
end;



Recursmat:=function(S,p,B)
local i, j, k, vect, eig, t, temp2, temp,SE;
temp:=[];
i:=1;
t:=0;
eig:=[];
k:=Size(S);
while i<=k do
	if S[i]=[] then
	Print("");
	else
	eig[i]:=Eigenvalues(GF(p),S[i]);
	if Size(eig[i])=Size(S[i]) then
	Print("SURPRISE \n ");
	t:=i;
	i:=k;
	fi;
	fi;
	i:=i+1;

od;
Print("t = ",t,"\n");
if t=0  then
	Print("We have to use recursion? \n");
	#None of the matrices S[i] had distinct eigenvalues so we 
	#run Eigensplit on this set S and then with the answer we make our vectors
	Print("What Shall We do here \n\n");
	SE:=EmptyRemove(S);
	Print("SE = ",SE,"\n\n\n");
	temp:=EigenSplit(SE,p);
	Print("temp = ",temp,"\n");
	#Temp is a set of liearly indepentent vectors, which we use to get 	our vectors
	if temp=[] then
		vect:=[];
	else
		Print("B = ",B,"\n");
		vect:=Recovereig(temp,p,B);
		Print("vect = ",vect,"\n");
	fi;
else
Print("We can get some vectors here \n ");
vect:=SplitMatrixEig(B, p, S[t]);
fi;

return vect;
end;

EmptyRemove:=function(M)
local i, L, k;
i:=1;
L:=[];
k:=Size(M);
while i<=k do
	if M[i] = [] then
	Print("");
	else
	Add(L,M[i]);
	fi;
	i:=i+1;
od;
return L;
end;

#This function has a bug and so doesn't wok for all M#
#Given a set of commuting matrices this will find the set of linearly independent eigenvectors which are eigenvectors for all the matrices.  Assume that the first matrix M[1] is the identity, p is the degree of the galois field we are working over.  #

EigenSplit:=function(M, p)   
local i, j, s, eig, eigval, eigvect,eigsp ,k ,eigs, r, b, N, Nu, V, Diag, match, l, temp, q, coeffs, x, lininvect, Md, kd , eig2, iter, jk, secondlevel, temp4;
eig2:=[];
eigval:=[];
eigvect:=[];
eigsp:=[];
Diag:=[];
N:=[];
i:=1;
k:=Size(M);
i:=1;
iter:=[];

#k is the number of Conjugcay classes - the number of matrices M

while i<=k do
	j:=1;
	eigval[i]:=Eigenvalues(GF(p),M[i]);
	eigvect[i]:=Eigenvectors(GF(p),TransposedMat(M[i]));	
	Diag[i]:=Inverse(TransposedMat(eigvect[i]))*M[i]*TransposedMat(eigvect[i]);
	Diag[i]:=ZToMat(Diag[i],p);
	#Print("Diag[",i,"] = ",Diag[i],"\n");
	i:=i+1;
	od;
i:=1;
while i<=k do
	eigval[i]:=ZToMat([eigval[i]], p);
	eigval[i]:=eigval[i][1];
	i:=i+1;
od;

s:=[];
i:=1;
while i<=k do
	s[i]:=Size(eigval[i]);
	i:=i+1;
od;

match:=EigvalEigvectMatch(eigval, eigvect, Diag, M, p, k, s);
iter:=[];

#We just choose the second matrix to bust open and then try the third etc until we have k eigenvectors#
i:=2;
lininvect:=[];
while i<=2 do
	Print("Number i = ",i,"\n");
	iter[i]:=[];
	jk:=1;
	while jk<=s[i] do
	iter[i][jk]:=[];
	Print("eigspace?:= ",jk,"\n");
	jk:=jk+1;
	od;
	V:=[];
	l:=1;
	while l<=k do;
		j:=1;	
		V[l]:=[];
		if l = i then
		Print("No point in splitting here! \n");
		fi;
		#Print("l = ",l,"\n");
		#l  loops over matrices M[l]#
		#Print("s[",i,"] = ",s[i],"\n");
		q:=0;
		while j<=s[i] do
			r:=1;
			
			V[l][j]:=[];
			#This loops over the basis of the eigenspaces of M[i]#
			#Print("match[",i,"][",j,"] = ",match[i][j]);
			#Print("Size match[",i,"][",j,"] =",Size(match[i][j]),"\n");
			while r<=Size(match[i][j]) do
				#Print("match[",i,"][",j,"][",r,"] = ",match[i][j][r]," \n");
				
				temp:=M[l]*match[i][j][r];
				temp:=ZToMat([temp],p);
				V[l][j][r]:=temp[1];
				r:=r+1;
			od;
			
			j:=j+1;
		od;
		l:=l+1;
	od;
	#V is the result of all these multiplications#
	Print("V = ",V,"\n");
	coeffs:=[];
	j:=1;
	q:=0;
	while j<=Size(V) do
		Print("j = ",j,"\n");
		l:=1;
		coeffs[j]:=[];
		eig2[j]:=[];
		while l<=Size(V[j]) do
			Print("match[",i,"][",l,"] = ",match[i][l],"\n");
			Print("V[",j,"][",l,"] = ",V[j][l]," \n");
			Print("l = ",l,"\n");
			if Size(lininvect) =k then
				l:=Size(V[j]);
				i:=k;
				Print("We have found all of our vectors \n");
				j:=Size(V);
			else
				temp:=SplitMatrix(V[j][l],match[i][l], p);
				if temp[1]=1 then
					coeffs[j][l]:=temp[2];
					Print("Success We can split at the moment \n");
				elif temp[1]=0 then
					if temp[2]=[] then
						Print("temp[2] is null");
					else 
						Print("l = ",l,"\n");
						Print("temp[2] =",temp[2],"\n");
					
						Add(iter[i][l],temp[2]);
					fi;
					Print("We are unable to split at the moment \n");
					coeffs[j][l]:=[];
				fi;
				r:=1;
				while r<=Size(coeffs[j][l]) do
					#Print("coeffs[",j,"][",l,"][",r,"] = ",coeffs[j][l][r],"\n");
					Add(lininvect,coeffs[j][l][r]); 
					r:=r+1;
				od;
			fi;
			l:=l+1;
			lininvect:=Set(lininvect);
		od;
	j:=j+1;
	od;
	#Print("iter[",i,"] = ",iter[i],"\n");#
i:=i+1;
od;
secondlevel:=[];
if Size(lininvect) =k  then
	Print("We have found all of our vectors \n"); 
else 
	#Print("iter[2][1][1] = ",iter[2][1][1],"\n");
	temp4:=Eigenvalues(GF(p),iter[2][1][1]);
	secondlevel[1]:=EigenSplit(iter[2][1], p, Size(iter[2][1]));
	#Print("Not all vectors found \n");
	Add(lininvect,secondlevel);
	
fi; 
	
lininvect:=ZToMat(lininvect,p);
lininvect:=Set(lininvect);
Print("lininvect = ",lininvect,"\n");
Print("coeffs =",coeffs,"\n");
Print("s = ",s,"\n");
return lininvect;

end;


Recovereig:=function(l,b,p)
local i, j, k, vect, temp;
k:=Size(l);
i:=1;
vect:=[];
while i<k do
	temp:=[];
	while j<=Size(l[i]) do
	temp:=temp+l[i][j]*b[j];
	j:=j+1;
	od;
	vect[i]:=ShallowCopy(temp);
i:=i+1;
od;
return vect;
end;

SplitMatrixEig:=function( b, p, split)
local i , j, k, c, q, eig, eigval ,eigvect, temp, r;
eigvect:=[];
eig:=Eigenvectors(GF(p),TransposedMat(split));
#Print("eig = ",ZToMat(eig,p),"\n");
#we now take the eigenvectors of this matrix and get the eigenvectors we are after#
eigvect:=[];
i:=1;
while i<=Size(eig) do
	j:=1;
	eigvect[i]:=[];
	temp:=0;
	while j<=Size(eig[i]) do
	r:=1;
		#Print("eig[",i,"][",j,"] =",Int(eig[i][j]),"\n");
		#Print("b[",j,"] = ",b[j],"\n");
		temp:=temp+b[j]*eig[i][j];
		j:=j+1;
	od;
	eigvect[i]:=temp;
	i:=i+1;
od;
#Print("eigvect:=",eigvect,"\n");	
return eigvect;
end;

FormSplitMatrix:=function(a,b,p)
local i , j, k, c, q, split, eig, eigval , eigvect, temp, r, check;
q:=[];
#Print("a =",a,"\n");
#Print("b =",b,"\n");
#Print("p =",p,"\n"); 
split:=[];
check:=0;
q[1]:=[];
q[2]:=[];
#we need to get the coefficients#
#a[i] is a multiple of some b[j]#
i:=1;
while i<=Size(a) do
	j:=1;
	while j<=Size(b) do
		k:=1;
		check:=0;
		while k<p do
			c:=b[j]*k*Z(p)^0;
			#Print("b[",j,"] = ",b[j],"\n");
			#Print("k =",k,"\n");
			#Print("c = ",ZToMat([c],p),"\n");
			#Print("a[",i,"] = ",a[i],"\n");
			if c = a[i]*Z(p)^0 then
				#Print("We have a match \n");
				check:=1;
				q[1][i]:=k;
				q[2][i]:=j;
				#Print("q = ",q,"\n");
				k:=p+3;
				j:=Size(b);
			else
				k:=k+1;
			fi;
		od;
	j:=j+1;
	od;
	if check=0 then
		Print("This matrix is rubbish \n ");
		i:=Size(a);
	fi;
i:=i+1;
od;
#Print("q = ",q,"\n");

if Size(q[1])=Size(a) and check=1 then
	i:=1;
	while i <=Size(q[1]) do
		j:=1;
		split[i]:=[];
		while j<=Size(q[1]) do
			split[i][j]:= 0;
		j:=j+1;
		od;
		i:=i+1;
	od;
	i:=1;
	while i<=Size(q[1]) do;
		split[q[2][i]][i]:=q[1][i];
		i:=i+1;
	od;
	split:=split*One(GF(p));
else
	Print("This eigenspace can't be broken up with this matrix \n");
	split:=[];
fi;
#Print("split =", ZToMat(split,p),"\n");
return split;
end;


SplitMatrix:=function(a, b, p)
local  split, eigvect, eigval;
eigvect:=[];
split:=FormSplitMatrix(a,b,p);
if split=[] then
	eigvect[1]:=0;
	eigvect[2]:=split;
else
eigval:=Eigenvalues(GF(p),split);
Print("Eigenvalues of split are ",eigval,"\n");
if Size(eigval) = Size(split) then
	eigvect[1]:=1;
	eigvect[2]:=SplitMatrixEig(b, p, split);
else
	Print("Split doesn't have distinct eigenvalues \n");
	eigvect[1]:=0;
	eigvect[2]:=split;
fi;
fi;
return eigvect;

end;			
 
EigvalEigvectMatch:=function(eigval, eigvect, Diag, M, p, k, s)
local i, j,l, match, q, z;
match:=[];
i:=1;
while i<=k do
	j:=1;
	#Print("Diag[",i,"] = ",Diag[i],"\n");
	match[i]:=[];
	while j<= s[i] do 
		l:=1;
		q:=0;
		match[i][j]:=[];
		while l<=k do
			#Print("Diag[",i,"][",l,"][",l,"] = ", Diag[i][l][l], "\n");
			#Print("eigval = ",eigval[i][j],"\n");
			if Diag[i][l][l]=eigval[i][j] then
			q:=q+1;
			#Print("eigvect[",i,"][",l,"] = ",ZToMat([eigvect[i][l]],p),"\n");
			z:=ZToMat([eigvect[i][l]],p);
			match[i][j][q]:=z[1];
			#Print("match[",i,"][",j,"][",q,"] = ",match[i][j][q],"\n");
			fi;
		l:=l+1;
		od;
	j:=j+1;
	od;
#Print("match[",i,"] = ",match[i],"\n");
i:=i+1;
od;

#Print("match = ",match,"\n");
return match;

end;



#This may find a set of simultaneous eigenvectors - though it is not very reliable for large p
EigvectSearch:=function(M,p)
local i , j, k, R, a, b, c, d, s;
R:=0;
k:=Size(M);
i:=1;
M:=M*One(GF(p));
j:=1;
while j <= 30 do
		R:=[];
		i:=1;
		while i <=k do
		c:=Random(0,p);
		Print("c:= ",c,"\n");
		R:= R + c*M[i];
		i:=i+1;
		od;
a:=Eigenvalues(GF(p),R);
b:=Eigenvectors(GF(p),TransposedMat(R));
d:=b*R*TransposedMat(b);
#Print("d = ",d,"\n");
#Print("a = ",a,"\n");
#Print("there are ",Size(a)," eigenvalues \n");
if Size(a) = k then
j:=301;
else
j:=j+1;
fi;
od;

return b;
end;			
 

#This is the function which returns the character table of a group, try not to put "big" groups in.  due to a bug in EigenSplit (a lack of enough recursion) This might not 
# work for your group.  If it does return an answer it is the correct one.

DixonAlg:=function(G)
local i, j, k, m, a, b, c, d, r, p, q, V, U, H,K, M, h, g, s, Diag, eig, herm, eigval, deg, char, z, pow, l, ans, eig2, charcom, P;
#K:=[];
if G=Group(()) then
	charcom:=[];
	else
K:=HConjClass(G,G);
H:=[];
d:=[];
deg:=[];
ans:=[];
pow:=[];
char:=[];

h:=[];
k:=Length(K);
#Print("k = ",k,"\n");
#Find the sizes of the K-conjugacy classes#
i:=1;
while i <=k do
	h[i]:=Size(K[i]);
	i:=i+1;
od;

#Finds a suitable prime (only up to 200000 at the moment)#
m:=Exponent(G);
#Print("Exponent = ",m,"\n");
g:=Size(G);
q:= 4*g;
i:=1;
while i<=168 do
	p:=Primes[i];
	if p*p >q then
		q:=p-1;
		if q mod m = 0 then
			#Print("p = ",p,"\n");
			i:=200;
		fi;
	fi;
	i:=i+1;
od;


if i=169 then
i:=997;
while i<=200000 do
	if IsPrimeInt(i) then
		#Print("the no ",i," is prime\n");
		p:=i;
		if p*p >q then
			q:=p-1;
			if q mod m = 0 then
			#	Print("p = ",p,"\n");
				i:=200001;
			fi;
		fi;
	fi;
	i:=i+1;
od;
fi;
i:=1;
M:=ClassCoeffs(K);

#Print("M:= ",M,"\n");

while i<=k do
	M[i]:=M[i]*One(GF(p));
	i:=i+1;
od;

eig:=[];
r:=1;
while r<=k do
	eig[r]:=Eigenvectors(GF(p),TransposedMat(M[r]));
r:=r+1;
od;


i:=2;
s:=0;
eigval:=[];
while i<=k do
	eigval[i]:= Eigenvalues(GF(p),M[i]);
i:=i+1;
od;


i:=2;
while i<=k do
	#Print("Sizeeigval[",i,"] =",Size(eigval[i]),"\n");
	if Size(eigval[i])= k then
		#Print("we'll use the eigenvectors from matrix M[",i,"] \n");
		#Print("eig[",i,"] = ",ZToMat(eig[i],p),"\n");#
		s:=i;
		eig2:=eig[s];
		i:=k+1;
	else
		i:=i+1;
	fi;
od;

if s=0 then
	#Print("Couldn't find any linearly independent eigenvectors \n"); 
	#Print("we will have to split the eigenspace of M[2] \n");\
	eig2:=EigenSplit(M, p);
	#eig2:=EigvectSearch(M,p);
		
	eig2:=eig2*One(GF(p));
	#Print("eig2 = ",eig2,"\n");
fi;

	i:=2;
	Diag:=[];
	while i<=k do
	#	Print("eig[",i,"] =",eig[i],"\n"); 
	#	Print("eig[",i,"] =",ZToMat(eig[i],p),"\n"); 
	#	Print("M[",i,"] = ",M[i],"\n");
		Diag[i]:=Inverse(TransposedMat(eig2))*M[i]*TransposedMat(eig2);
	#	Print("Diag[",i,"] = ",Diag[i]," \n");
	#	Print("Diag[",i,"] = ",ZToMat(Diag[i],p)," \n");
	#	Print("Rank eig[",i,"] = ",Rank(TransposedMat(eig[i])),"\n");
	#	Print("eigval[",i,"] = ",eigval[i],"\n");
	#	Print("eigval[",i,"] = ",ZToMat([eigval[i]],p),"\n");
		i:=i+1;
	od;
deg:=ComputeDegDixon(eig2, p, G, K);
char:=ComputeCharDixZp(deg,eig2, p, G, K);
charcom:=ComputeCharComplex(char, deg, p, G, K, m);
fi;
return charcom;
end;

#This function sets the trivial character as \chi_1#
OneChar:=function(c, k)
local i, j, char, trivchar, q, temp;
trivchar:=[];
i:=1;
while i<=k do
	trivchar[i]:=1;
	i:=i+1;
	od;
char:=[];
if c[1]=trivchar then
	#Print("Triv char already in first place \n");#
	q:=1;
else
	q:=0;
fi;
i:=1;
if q=0 then
	while i<=k do	
		if c[i]=trivchar then
		#Print("trivchar in the ",i," place \n");#
		temp:=c[1];
		c[1]:=c[i];
		c[i]:=temp;
		i:=k;
		fi;
	i:=i+1;
	od;
fi;


end;

#Find a suitable z#
		
Computez:=function(char, deg, p, G, H, m)
local i, j, k, sum, d, a, ci, dzqzi, q, dz, qz, qzi, r, rp, ch, prim, z, ans, iz, ex, f, pow, charac, n, s, w;
k:=Size(H);
ch:=[];
z:=0;
i:=1;
q:=1;
#Here we find a suitable integer z#
while i<=p do;
	j:=i*Z(p)^0;
	if Order(j)= m then
		z:=j;
		i:=p+1;
	else
		i:=i+1;
	fi;
od;

#Print("z=",z,"\n");
#Print("Int(z) = ",Int(z),"\n");

return z;
end;

#Given the characters in Z_p  we now find the characters in C#

ComputeCharComplex:=function(char, deg, p, G, K, m)
local i, j, k, sum, d, a, ci, dzqzi, q, dz, qz, qzi, r, rp, ch, prim,  z, ans, iz, ex, f, pow, charac, n, s, w;
k:=Size(K);
ch:=[];
q:=1;
prim:=E(m);
z:=Computez(char, deg, p, G, K, m);
#Print("z=",z,"\n");
#Print("Int(z) = ",Int(z),"\n");

n:=0;
pow:=[];
while n<m do;
	pow[n+1]:=[];
	pow[n+1]:=PowerConjClassDash(K,n);
	n:=n+1;
od;
#Print("pow:= ",pow,"\n");	
i:=1;
j:=1;
s:=0;
n:=1;
sum:=0;
w:=[];
while i<=k do;
	j:=1;
	w[i]:=[];
	while j<=k do;
		s:=0;
		w[i][j]:=[];
		while s<m do;
			n:=0;
			sum:=0;
			while n < m do;
				sum:=sum+char[i][pow[n+1][j]]*(z^(-s*n));
				n:=n+1;
			od;
			w[i][j][s+1]:=sum*(1/m);
			w[i][j][s+1]:=Int(w[i][j][s+1]);
			s:=s+1;
		od;
	j:=j+1;
	od;
i:=i+1;
od;

#Print("w = ",w,"\n");

i:=1;
j:=1;
charac:=[];
sum:=0;
while i<=k do;
	j:=1;
	charac[i]:=[];
	while j<=k do;
		s:=0;
		r:=Order(K[j][1]);
		sum:=0;
		while s<m do;
			sum:=sum+w[i][j][s+1]*prim^s;
			s:=s+1;
		od;
		charac[i][j]:=sum;
		j:=j+1;
	od;
i:=i+1;
od;

#Make sure that \chi_1 is the trivial character#
OneChar(charac, k);

#Print("charac =", charac, "\n");
			
return charac;
end;

#Compute the characters in Z_p#

ComputeCharDixZp:=function(deg,eig, p, G, K)
local i, j, k, sum, d, ci, dzqzi, q, dz, qz, qzi, r, rp, char;
i:=1;
j:=1;
k:=Size(K);
char:=[];
while j<=k do;
	i:=1;
	char[j]:=[];
	while i<=k do;
	char[j][i]:=eig[j][i]*deg[j]*Inverse(Size(K[i])*Z(p)^0);
	char[j][i]:=Int(char[j][i]);
	i:=i+1;
	od;
j:=j+1;
od;

#Print("char =",char,"\n");

return char;
end;


#Given some linearly independent eigenvectors we compute the degrees of the representations in Z_p#
ComputeDegDixon:=function(eig, p, G, K)
local i, j, k, sum, deg, d, ci, dzqzi, q, dz, qz, qzi, r, rp, jd;
i:=1;
jd:=InvConjjdash(K);
#Print("jd = ",jd,"\n");
sum:=[];
r:=[];
rp:=[];
while i<=Size(eig) do;
	#Print(" i = ",i,"\n");
	sum[i]:=HermitInnerProd(ZToMat([eig[i]],p),ZToMat([eig[i]],p),G,K,p,jd);
	#Print("ZToMat(eig[",i,"]) = ",ZToMat([eig[i]],p),"\n");#
	#Print("sum[",i,"] = ",sum[i],"\n");#
	i:=i+1;
od;

Print("sum =",sum,"\n");
i:=1;
while i<=Size(sum) do;
	r[i]:=[];
	d:=DenominatorRat(sum[i]);
	q:=sum[i]*d;
	dz:= d*Z(p)^0;
	qz:= q*Z(p)^0;
	qzi:=Inverse(qz);
	dzqzi:=dz*qzi;
	ci:=Int(dzqzi);
	r[i]:=RootsMod(ci,p);
	Print("r[",i,"] = ",r[i],"\n");
	if r[i][1]  < p/2 then
	rp[i]:=r[i][1];
	else
	rp[i]:=r[i][2];
	fi;
i:=i+1;
od;
	
	
#Print("r =",r,"\n");	
Print("deg =",rp,"\n");
return rp;
end;


#Given a conjugacy class this returns the inverse of this conjugacy class#
InverseConjclass:=function(K)
local i, j, I;
i:=1;
I:=[];
j:=1;
while j<=Size(K) do;
	i:=1;
	I[j]:=[];
	while i<=Size(K[j]) do;
	I[j][i]:=Inverse(K[j][i]);
	i:=i+1;
	od;
j:=j+1;
od;
return I;
end;

#Computes a Hermitian inner product described by a lemma of Burnsides#
HermitInnerProd:=function(u,w,G,K,p,jd)
local i, j, k, sum;
sum:=Z(p)*0;
#Print("u =",u,"\n");
#Print("w =",w,"\n");
i:=1;
j:=1;
while i<=Size(u[1]) do;
		sum:= sum+((u[1][i]*w[1][jd[i]]*Z(p)^0)*Inverse(Size(K[i])*Z(p)^0));
	i:=i+1;
od;
sum:=sum*Inverse(Size(G)*Z(p)^0);
sum:=Int(sum);

return sum;
end;

#THis uses Burnsides's method to compute the character values but only if there is a class coefficient matrix whose distinct
 #eigenvalues can be computed over the rationals and then the linearly independent eigenvectors can be computed over the rationals,
#We haven't done eigensplit for Burnside
Burnside:=function(G)
local i, j, k, x, y, z,l, m, n, h,  H, r, s, t, c, p , q , g, C, eigval, M, e, pq, eig, V, B, d, D, Grp, Diag, sp;
M:=[];
H:=HConjClass(G,G);
M:=ClassCoeffs(H);
h:=[];
H:=[];
c:=[];
D:=[];

eigval:=[];
k:=Length(H);
Print("k = ",k, "\n");
#Find the sizes of the H-conjugacy classes#
i:=1;
while i <=k do
	h[i]:=Size(H[i]);
	i:=i+1;
od;

p:=1;
q:=1;
r:=1;
eig:=[];
r:=1;
while r<=k do
	eig[r]:=Eigenvectors(Rationals,TransposedMat(M[r]));
	r:=r+1;
od;
e:=[];

i:=2;
s:=0;
eigval:=[];
while i<=k do
eigval[i]:= Eigenvalues(Rationals,M[i]);
i:=i+1;
od;
i:=2;

while i<=k do
	Print("Sizeeigval[",i,"] =",Size(eigval[i]),"\n");
	if Size(eigval[i])= k then
		Print("we'll use the eigenvectors form matrix M[",i,"] \n");
		Print("eig[",i,"] = ",eig[i],"\n");
		s:=i;
		i:=k+1;
	else
		i:=i+1;
	fi;
od;

if s=0 then
	Print("Couldn't find any linearly independent eigenvectors \n"); 
else
	i:=2;
	Diag:=[];
	while i<=k do
		Print("eig[",i,"] =",eig[i],"\n"); 
		Diag[i]:=Inverse(TransposedMat(eig[s]))*M[i]*TransposedMat(eig[s]);
		Print("Diag[",i,"] = ",Diag[i]," \n");
		Print("Rank eig[",i,"] = ",Rank(TransposedMat(eig[i])),"\n");
		Print("eigval[",i,"] = ",eigval[i],"\n");
		i:=i+1;
	od;
	D:=ComputeDegs(eig[s],G,H);
fi;



return D;
end;

# Given a group and a set of eigenvectors B which spans the eiegenspace of the 
# Matrices M_r we can calculate the degrees of the irreducible characters#

ComputeDegs:=function( W, G, K)
local i, j, r, u, v, k, CharDeg, c, Char, rdash;
k:=Size(K);
rdash:=[];
rdash:=InvConjjdash(K);
CharDeg:=[];
Char:=[];
i:=1;
while i<=k do
	r:=1;
CharDeg[i]:=0;	
	while r<=k do
		CharDeg[i]:=CharDeg[i]+ W[i][r]*W[i][rdash[r]]/(Size(K[r]));
		r:=r+1;
	od;
i:=i+1;
od;

i:=1;
while i<=k do
	CharDeg[i]:=RootInt(Size(G)/CharDeg[i]);
	i:=i+1;
od;
i:=1;
j:=1;
Char:=[];
while i<=k do
	j:=1;
	Char[i]:=[];
	while j<=k do;
		Char[i][j]:=W[i][j]*CharDeg[i]/(Size(K[j]));
		j:=j+1;
	od;
i:=i+1;
od;

#Print("Char =",Char,"\n");

#Print("CharDeg =",CharDeg,"\n");
return Char;
end;


SortedConjClass:=function(G,tbl)
local i, j, k, K, Conj1, Conj2, ident;
Conj1:=ConjugacyClasses(G);
k:=Size(Conj1);
Conj2:=[];
K:=[];
i:=1;
while i <= k do
Conj2[i]:=Elements(Conj1[i]);
i:=i+1;
od;
ident:=IdentificationOfConjugacyClasses(tbl);
#Print("ident:= ",ident,"\n");
K:=[];
i:=1;
while i<=k do
K[i]:=ShallowCopy(Conj2[ident[i]]);
i:=i+1;
od;


return K;

end;

ClassEqual:= function(L,P)
local i, j, k, p, q;
i:=1;
#Print("Size of L =",Size(L),"\n");#
#Print("Size of P =",Size(P),"\n");#
while i<= Length(L) do
	j:=1;
	#Print("Size of [",i,"] = ",Size(L[i]), "\n");#
	while j<= Length(P) do
		if IsEqualSet(L[i],P[j]) then
			j:= Length(P)+1;
			p:=1;
		else
			j:=j+1;
			p:=0;
		fi;
	od;
	if p=1 then
		i:=i+1;
	else
		q:=i;
		i:=Length(L)+1;
	fi;
od;

if p=1 then
	#Print("equal \n");#
	k:=1;
else
	#Print("not equal in set ",q,"\n");#
	k:=0;
fi;

return k;

end;




DixonAlg1:=function(G)
local tbl,char;
tbl:=CharacterTable(G);;
char:=List(Irr(tbl),ValuesOfClassFunction);;
return char;
end;

DixonAlg3:=function(tbl)
local char;
char:=List(Irr(tbl),ValuesOfClassFunction);;
return char;
end;

DixonAlg2:=function(G)
local i, j, k, char, tbl;
tbl:=CharacterTable(G);;
return tbl;
end;