LoadPackage("chop");
Read("/u6/jramalheiratsu/SchurCoverOfSymmetricGroupConjugacyClasses");

#Finds the lcm of a list.
LcmList := function(list)
local l, lcm, i,j, counter;
l := Collected(list);
#check if all things in list are 1
if Length(l) = 1 then
	return 1;
fi;
#don't need to take lcm of (1,n)
if l[1][1] = 1 then
	Remove(l,1);
fi;
#keep taking pairwise lcm until you get one number
counter := Length(l);
while counter > 1 do
	lcm := [];
	for i in [1..Length(l)-1] do
		for j in [i+1..Length(l)] do
			Add(lcm, LcmInt(l[i][1],l[j][1]));
		od;
	od;
	l := Collected(lcm);
	counter := Length(l);
od;
return l[1][1];
end;

CMatz := function(mat,field)
local p,d,cf;
cf := Collected(Factors(Size(field)));
p := cf[1][1];
d := cf[1][2];
return CMat(List(mat, x -> CVec(x,p,d)));
end;

changefield := function(cgens,field)
local i,mats, p, basefield, n, m, lcm;
mats := [];
basefield := Field(cgens[1][1]);
p := Characteristic(field);
if not Characteristic(basefield) = p then
	Error("cannot change field as the specified field has a different characteristic");
fi;
n := Dimension(basefield);
m := Dimension(field);
lcm := LcmInt(m,n);
for i in [1..Size(cgens)] do
	Add(mats,CMat(List(cgens[i],x -> CVec(Unpack(x),GF(p^lcm)))));
od;
return mats;
end;

#Input: Group, prime, list of chars (their number in GAP).
#Output: the minimum field all chars are over

minfield := function(G,p,l)
local list,fields,ct, irr, i,j;
list := [];
fields := [];
ct := CharacterTable(G) mod p;
irr := Irr(ct);
for i in l do
	Add(list, List(irr[i],x->FrobeniusCharacterValue(x,p)));
od;
for j in [1..Length(l)] do
	Add(fields, LogInt(Size(Field(list[j])),p));
od;
return GF(p^LcmList(fields));
end;

#Input: Two lists of matrices. kprod takes the tensor product of the lists of matrices. Output: A list of matrices. (Note matrices should be as CMat).
#add a check that matrices are over the same field.

kprod := function(s,t)
local i,kp;
kp := [];
for i in [1..Length(s)] do
	Add(kp,KroneckerProduct(s[i],t[i]));
od;
return kp;
end;

#-------------------------------------------------------------

trivialrep := function(n)
local i, tr;
tr := [];
for i in [1..n] do
	Add(tr,1);
od;
return tr;
end;

mosttensor := function(ct,tensorlist)
local irr, i,j,k, counter, tensors, sortten, mult, max, maxten,degree, smallest, remaining;
irr := Irr(ct);
smallest := [];
counter := Filtered(tensorlist, x -> Length(x)>0);
tensors := [];
k := 1;
while Length(counter) > 0 do
	remaining := [];
	tensors := [];
	for i in [1..Length(counter)] do
		for j in [1..Length(counter[i])] do
			Add(tensors,counter[i][j]);
		od;
	od;
	sortten := Collected(tensors);
	mult := List(sortten,x->x[2]);
	max := Positions(mult,Maximum(mult));
	maxten := sortten{max}[1];
	degree := List(maxten, x-> irr[x[1]][1]*irr[x[2]][1]);
	Add(smallest, maxten[Position(degree,Minimum(degree))]);
	for i in [1..Length(counter)] do
		if not smallest[k] in counter[i] then
			Add(remaining,i);
		fi;
	od;
	counter := counter{remaining}; 
	k := k+1;
od;
return smallest;
end;

generalplan := function(ct,knownchars)
local irr, tr, step, i, factors, kn, j,k, allreps, tensor, decomp, newchars, tensorlist, plan, notmade;

kn := [];
for i in [1..Length(knownchars)] do
	Add(kn,knownchars[i]);
od;
irr := Irr(ct);
tr := Position(irr,trivialrep(Size(irr)));

#check to make sure knownchars is not trivial
if tr in kn then
	Remove(kn,Position(kn,tr));
fi;
if Length(kn) = 0 then
	Error("Must give a nontrivial representation.");
fi;

#create list of which representations you need,remove trivial
factors := [];
for i in [1..Size(irr)] do
	if not i in kn then
		Add(factors,i);
	fi;
od;

allreps := false;
newchars := kn;
plan := [];
while allreps = false do
	#clear out tensorlist; use this to keep track of which 	#tensors contain each irreducible
	tensorlist := [];
	for i in [1..Size(irr)] do
		Add(tensorlist,[]);
	od;

	#find which irreducibles occur in each tensor
	for i in kn do
		for j in newchars do
			tensor := Tensored([irr[i]],[irr[j]]);
			decomp := Decomposition(irr,tensor,5);
			for k in factors do
				if decomp[1][k]>0 then
					Add(tensorlist[k],[i,j]);
				fi;
			od;
		od;
	od;

	#figure out smallest set of tensors
	Append(plan,mosttensor(ct,tensorlist));

	#record the new characters you have,remove them from factors
	newchars := [];
	for k in [1..Length(tensorlist)] do
		if Length(tensorlist[k]) > 0 then
			Add(newchars,k);
			Remove(factors,Position(factors,k));
		fi;
	od;
	Append(kn,newchars);
	if Length(newchars) = 0  or Length(kn) = Length(irr)-1 then
		allreps := true;
	fi;
od;

#check if all characters were created
if Length(kn) < Length(irr)-1 then
	notmade := [];
	for i in [1..Length(irr)] do
		if not i in kn and not i = tr then
			Add(notmade,i);
		fi;
	od;
	Print("The following characters cannot be constructed: ");
	Print(notmade);
fi;

return plan;
end;

SchurVal := function(n,p,cgens)
local irr, classes, pclasses, traces;
#convert char table to Jon's char table
#irr := Irr(CharacterTable(G) mod p);
irr := Irr(CharacterTable("DoubleCoverSymmetric",n) mod p);
classes := SchurCoverOfSymmetricGroupConjugacyClasses(n,cgens);
pclasses := Filtered(classes, x->Order(x) mod p>0);
traces := List(pclasses, x->BrauerCharacterValue(Unpack(x)));
return traces;
end;

SchurChar := function(n,p,cgens)
local irr, classes, pclasses, traces, tr, dims;
irr := Irr(CharacterTable("DoubleCoverSymmetric",n) mod p);
classes := SchurCoverOfSymmetricGroupConjugacyClasses(n,cgens);
pclasses := Filtered(classes, x->Order(x) mod p>0);
traces := List(pclasses, x->BrauerCharacterValue(Unpack(x)));

#Deal with the fact that 1 dim chars will not remove the pclasses
if traces[1] = 1 then
	#check to make sure that there are only two 1 dim reps
	dims := irr{[1..Size(irr)]}[1];
	if Length(Positions(dims,1))>2 then
		Error("Problem with SchurChar on the 1 dim chars");
	fi;
	#Figure out which 1 dim rep it is
	tr := trivialrep(Length(traces));
	if traces = tr then
		return 1;
	else
		return Positions(dims,1)[2];
	fi;
fi;

return Position(irr,traces);
end;

SchurSplittingField := function(n,p)
local list,powers, irr, i,j;
list := [];
powers := [];
irr := Irr(CharacterTable("DoubleCoverSymmetric",n) mod p);
for i in [1..Length(irr)] do
	Add(list, List(irr[i],x->FrobeniusCharacterValue(x,p)));
od;
for j in [1..Length(list)] do
	Add(powers, LogInt(Size(Field(list[j])),p));
od;
return GF(p^LcmList(powers));
end;


SchurDecomposeTensorProd := function(n,p,cgens1,cgens2)
local field, tensor, factors, irreps, chars, rep1,rep2,reps;
field := SchurSplittingField(n,p);
rep1 := changefield(cgens1,field);
rep2 := changefield(cgens2,field);
tensor := Module(kprod(rep1,rep2));
factors := Chop(tensor);
reps := List(factors.db, x-> List(RepresentingMatrices(x), y -> CMat(y)));
chars := List(reps, x-> SchurChar(n,p,x));
return rec(representations := reps, characters := chars);
end;

SchurFindAllReps := function(n,p)
local ct, irr, gpname, dims, cgens, reps, chars, plan, cgens1, cgens2, decomp, reprec, i, j;

ct := CharacterTable("DoubleCoverSymmetric",n) mod p;
irr := Irr(ct);
dims := irr{[1..Size(irr)]}[1];
gpname := Concatenation(["2.S",String(n)]);

#get starting representation and character
cgens := BasicSpinRepresentationOfSymmetricGroup(n,p,+1);
cgens := List(cgens, x-> CMat(x));
Add(cgens,cgens[1]^2);
reps := [cgens];
chars := [SchurChar(n,p,cgens)];

#get the plan
plan := generalplan(ct,chars);

#implement the plan
for i in [1..Length(plan)] do
	if Position(chars,plan[i][1]) = fail or
		Position(chars,plan[i][2]) = fail then
		Error("Something is wrong. You don't have the character that the plan wants you to use.");
	fi;
	cgens1 := reps[Position(chars,plan[i][1])];
	cgens2 := reps[Position(chars,plan[i][2])];
	decomp := SchurDecomposeTensorProd(n,p,cgens1,cgens2);
	
	#Add the new characters and representations to list
	for j in [1..Length(decomp.characters)] do
		if Position(chars,decomp.characters[j]) = fail then
			Add(reps, decomp.representations[j]);
			Add(chars, decomp.characters[j]);
		fi;
	od;
od;
#Now we should have all the representations, so just create the output record.
reprec := [];
for i in [1..Length(chars)] do
	Add(reprec, rec(
		character := chars[i],
		dim := dims[chars[i]],
		groupname := gpname,
		ring := Field(reps[chars[i]][1][1]),
		generators := reps[chars[i]]));
od;
return reprec;
end;