AttachSpec("/usr/local/magma/package/spec");

// The generator ordering in D5 is different in magma and in Fokko's program.
// This only the case only in D5, which I find very strange!
// For this reason we have to explicitly give the coxeter matrix for D5.

// MD5 := Matrix([[1,2,3,2,2],[2,1,3,2,2],[3,3,1,3,2],[2,2,3,1,3],[2,2,2,3,1]]);

// Change group type here:

// W<s,t,u> := CoxeterGroup(GrpFPCox, "H3");

// W<s,t,u,v,w> := CoxeterGroup(GrpFPCox, MD5);

W<s,t,u,v> := CoxeterGroup(GrpFPCox, "H4");

// W<s,t,u, v> := CoxeterGroup(GrpFPCox, "F4");

// W<s,t,u,v,w> := CoxeterGroup(GrpFPCox, "B5");

// W<s,t,u,v,w,x> := CoxeterGroup(GrpFPCox, "B6");

// W<s,t,u,v,w,x,y> := CoxeterGroup(GrpFPCox, "A7");

load "H4.m";

n := Rank(W);

CoxElt := [ W!CoxEltForms[i] : i in [1..#CoxEltForms]];

w0 := CoxElt[#CoxElt];

// There is one small error in the W-graphs that Fokko supplies:
// The indexes in the W-graph are one to small (because in a sequence 0 is the smallest for him, not 1).
// We fix this here:
// NOTE: Be careful not to reload this file without reloading the W-graph, 
// as calling this twice destroys all the numbering!
for i in [1..#WGraph] do
	for j in [1..#WGraph[i][2]] do
		WGraph[i][2][j][1] +:= 1;
	end for;
end for;

// To got back from the Weyl group to Fokko's indexing:
I := function(w)
	return Index(CoxElt, w);
end function;

// Kinder reports what elements break off in degree 0 when translate w down by s.
// (Note that we have to do everything upside down!)
Kinder := function(lr, m, s)
	i := m;
	d := s;
	
	if (lr eq "l" and d notin WGraph[i][1][1]) or (lr eq "r" and d notin WGraph[i][1][2])  then
		return {};
	end if;
	
	if lr eq "l" then
		return { WGraph[i][2][l][1]  : l in [1..#WGraph[i][2]] | d notin WGraph[WGraph[i][2][l][1]][1][1] };
	elif lr eq "r" then
		return { WGraph[i][2][l][1]  : l in [1..#WGraph[i][2]] | d notin WGraph[WGraph[i][2][l][1]][1][2] };
	end if;
	
	print "error: Kinder should be called with either \"l\" (left) or \"r\" (right) as lr paramater";
	return {};
end function;

// This function returns the elements which don't lie in \sigma(W) as well as the value of f_W on them.
// lower and upper may be set so as to only investigate part of the W-graph.

// The program has one major curiosity:
// We need to calculate constantly which elements should "break off" when we translate.
// Hence we sould expect and arrow x --> w in the W-graph when the KL basis element _H_w
// occurs in a product _H_s _H_w.
// Unfortunately, Fokko give us exactly the opposite - There is an arrow x --> w in the W-graph
// precisely when _H_x occurs in some product _H_s _H_w or _H_w _H_s.
// Hence we do everything ''upside down".
// Thus at the end we need to multiply all the elements by the longest element to get the real sets
// \sigma(W) and the critical pairs.

ProblemKinder := function(lower, upper)
	
	// A counter so that we have some idea how it is going:
	Z := ResidueClassRing(1000);
	counter := Z!0;
	
	// Problems records the elements which don't lie in \sigma(W).
	Problems := [];
	// Clumping records the values of the function f_W on them.
	Clumping := [];
	
	// We will go through using Fokko's ordering of the Coxeter elements.
	for i in [lower..Minimum(upper, #CoxEltForms)-1] do
		
		counter +:= 1;
		if counter eq Z!0 then
			printf " ... %o ... ", i;
		end if;
		
		w := CoxElt[i];
		Offspring := {};
		
		LeftAscent := [ d : d in [1..n] | d notin WGraph[i][1][1] ];
		RightAscent := [ d : d in [1..n] | d notin WGraph[i][1][2] ];
		
		for d in LeftAscent do
			
			j := Index(Problems, I((W.d)*w));
			
			// In the inductive definition of f_W we need to make sure that 
			// there exists at least one element sw or ws less than w such that
			// all elements in f_W(sw) or f_W(ws) are raised by s.
			// In fact we check something weaker: we check that there exists at
			// least one sw or ws such that f_W(sw) = {sw} (and similarly on the right).
			// If this is not the case we call it a "degreeZeroProblem":
			degreeZeroProblem := true;
			
			// what are the offspring?
			if Offspring eq {} then
				
				if j eq 0 then
					degreeZeroProblem := false;
					Offspring := Kinder("l", I((W.d)*w), d);
				else
					for Mod in Clumping[j] do
						Offspring join:= Kinder("l", Mod, d);
					end for;
				end if;
			else
				
				if j eq 0 then
					degreeZeroProblem:= false;
					Offspring meet:= Kinder("l", I((W.d)*w), d);
				else
					BigClump := {};
					for Mod in Clumping[j] do
						BigClump join:= Kinder("l", Mod, d);
					end for;
					Offspring meet:= BigClump;
				end if;
			end if;
			
			if (((Offspring eq {}) or (Offspring eq {i}))) and (degreeZeroProblem eq false) then
				continue i;
			end if;
		end for;
		
		for d in RightAscent do // right descents
			
			j := Index(Problems, I((W.d)*w));
		
			if j eq 0 then
				degreeZeroProblem := false;
				Offspring meet:= Kinder("r", I(w*(W.d)), d);
			else
				BigClump := {};
				for Mod in Clumping[j] do
					BigClump join:= Kinder("r", Mod, d);
				end for;
				Offspring meet:= BigClump;
			end if;
		
			if (Offspring eq {}) or (Offspring eq {i}) then
				printf "%o OK (but needed right). ", i;
				continue i;
			end if;
		end for;
		
		if degreeZeroProblem eq true then
			printf "POSSIBLE DEGREE ZERO PROBLEM!\n";
		end if;
		
		printf "\n%o: PROBLEM!!\n", i;
		printf "Problems: %o ... adding to sequence ... \n\n", Offspring;
		
		Append(~Problems, i);
		Append(~Clumping, Setseq(Offspring));
		
	end for;
	
	printf "finishing up ...\n\n";
	
	return Problems, Clumping;
	
end function;

// The right and left descents of elements of a Weyl group:
DescentL := function(w)
	return { W.i : i in [1..Rank(W)] | Length(W.i*w) le Length(w) };
end function;

DescentR := function(w)
	return { W.i : i in [1..Rank(W)] | Length(w*W.i) le Length(w) };
end function;


// w should be an element of the symmetric group.
// This routine returns the corresponding permutation.
permutation := function(w)
	a := Eltseq(w);
	perm := [];
	for i in [1..(n+1)] do
		perm[i] := i;
		for k in [1..#a] do
			if perm[i] eq a[k] + 1 then
				perm[i] -:= 1;
			elif perm[i] eq a[k] then
				perm[i] +:= 1;
			end if;
		end for;
	end for;
	return perm;
end function;

// Produces files of problem pairs:
FileOutput := procedure(a,b)
	printf "// Critical indices and pairs for the Coxeter group:\n";
	print W;
	printf "\n\n// %o of the %o elements of W do not lie in \\sigma(W).\n// They are:\n\n", #a, #W;

	printf "c := [ \n";
	for i in [1..(#a-1)] do
		printf "%o,\n", CoxElt[a[i]]*w0;
	end for;
	printf "%o \n ];", CoxElt[a[#a]]*w0;
	
	printf "\n\n\n // Values of f_W on the above elements (critical sets):\n\n p := [ \n";
		
	for i in [1..#b] do
		crit := a[i];
		printf "[ %o ", CoxElt[a[i]]*w0;
		for j in b[i] do
			if j ne crit then
				printf ", %o ", CoxElt[j]*w0;
			end if;
		end for;
		printf " ],\n";
	end for;

	printf " ];";

end procedure;

LatexOutput := procedure(ProblemString)
	
	AlreadyPrinted := {};

	for p in ProblemString do
		w1 := CoxElt[p[1][1]]*w0;
		w2 := CoxElt[p[2][1]]*w0;
		if Length(w2) ge Length(w1) then
			w3 := w1;
			w1 := w2;
			w2 := w3;
		end if;
		if w1^(-1) in AlreadyPrinted then
			continue p;
		end if;
		AlreadyPrinted join:= {w1};
		
		printf "$";
		for i in Eltseq(w1) do
			printf "s_%o", i;
		end for;
		printf "$ \n";
		printf "& $";
		for i in Eltseq(w2) do
			printf "s_%o", i;
		end for;
		printf "$ \\\\ \n";
	end for;
end procedure;

// For type A
LatexOutputA := procedure(ProblemString)
	
	AlreadyPrinted := {};

	for p in ProblemString do
		w1 := CoxElt[p[1][1]]*w0;
		w2 := CoxElt[p[2][1]]*w0;
		if Length(w2) ge Length(w1) then
			w3 := w1;
			w1 := w2;
			w2 := w3;
		end if;
		if w1^(-1) in AlreadyPrinted then
			continue p;
		end if;
		AlreadyPrinted join:= {w1};
		
		printf "$";
		for i in Eltseq(w1) do
			printf "s_%o", i;
		end for;
		printf "$ \n";
		printf "& $";
		for i in Eltseq(w2) do
			printf "s_%o", i;
		end for;
		printf "$ \\\\ \n";
		
	end for;
end procedure;

LatexOutputPermutation := procedure(ProblemString)
	
	AlreadyPrinted := {};
	
	for p in ProblemString do
		w1 := CoxElt[p[1][1]]*w0;
		w2 := CoxElt[p[2][1]]*w0;
		if Length(w2) ge Length(w1) then
			w3 := w1;
			w1 := w2;
			w2 := w3;
		end if;
		
		if w1^(-1) in AlreadyPrinted then
			continue p;
		end if;
		AlreadyPrinted join:= {w1};
		
				printf "$";
		for i in Eltseq(permutation(w1)) do
			printf "%o", i;
		end for;
		printf "$ ";
		printf "& $";
		for i in Eltseq(permutation(w2)) do
			printf "%o", i;
		end for;
		printf "$ \\\\ \n";
	end for;
end procedure;