// 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, "A3");

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

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

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

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

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

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

load "A7.m";

n := Rank(W);

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

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

// We now transform Fokko's data into something more suitable for us:
// Our W-Graph at w indexed by i is labelled by Vertices[i], and has edges going to Edges[i].
// Note: Because Fokko is working in C the W-graphs that Fokko supplies indexed 
// are one to small (because in a sequence 0 is the smallest for him, not 1).

Vertices := [];
Edges := [];

for i in [1..#WGraph] do
	Vertices[i] := [];
	Vertices[i][1] := WGraph[i][1];
	Edges[i] := [];
end for;

for i in [1..#WGraph] do
	Arrows := WGraph[i][2];
	for arr in Arrows do
		Edges[arr[1]+1] cat:= [i];
	end for;
end for;

// Kinder reports what elements break off in degree 0 when we translate w (i) up by s (d).
// Throughout we use the convention 1 = left, 2 = right.
Kinder := function(lr, i, d)
	
	// if d is in the left / right descent sets, then we just get (v + v^{-1}) * i.
	// will delete this check when I have more confidence in the program.
	if Index(Vertices[i][1][lr], d) ne 0 then
		printf "Index = %o, lr = %o: Shouldn't be here!", i, lr;
		return {};
	end if;
	return { l : l in Edges[i] | d in Vertices[l][1][lr] };	

end function;

// We now come to a series of routines designed to calculate \sigma(W) as well as the value of f_W on them.

// Those elements which do not belong to \sigma(W) (sorry for the name!).
Sigma := [];
// Clumping records the value of f_W on them.
Clumping := [];
// Blackhole records the elements in W where f_W is undefined.
Blackhole := [];

// returns the value of f_W at the index i, under the assumption that f_W is known on
// all elements less than i in the Bruhat order have been checked 
// (and therefore the sequences Sigma and Clumping accurately describe f_W).
// Note: Sigma and Clumping are passed to the routine, because magma was misbehaving
// when they were not past.
ClumpingAtIndex := function(i, Sigma, Clumping)

	// In the definition of \sigma(W) we need to find an element sw less than
	// w such that all elements of f_W(sw) are raised by sw. If this is not the
	// case then we have a degree zero problem.
	DegreeZeroProblem := true;
	
	// the guy we are dealing with
	w := W!(CoxEltForms[i]);
	
	IsFirstLoop := true;
	// Remember: 1 = left, 2 = right
	for lr in [1,2] do
		for d in Vertices[i][1][lr] do
			if lr eq 1 then
				j := I((W.d)*w);
			else
				j := I(w*(W.d));
			end if;
			
			l := Index(Sigma, j);
			if l eq 0 then
				source := [j];
				DegreeZeroProblem := false;
			else
				source := Clumping[l];
				// have to check that all elements in source are raised.
				// if this isn't the case we have to ignore this particular element.
				for jc in source do
					if Index(Vertices[jc][1][lr], d) ne 0 then
						continue d;
					end if;
				end for;
				DegreeZeroProblem := false;
			end if;
			Clump := {};
			for jc in source do
				Clump join:= Kinder(lr, jc, d);
			end for;
			
			if IsFirstLoop eq true then
				Offspring := Clump;
				IsFirstLoop := false;
			else
				Offspring meet:= Clump;
			end if;

		end for;
	end for;
	
	// output:
	if DegreeZeroProblem then
		printf "f_W not defined at %o\n", i;
		return {};
	else
		return Offspring;
	end if;
	
end function;

ProblemKinder := procedure(Anfang, Ende, ~Sigma, ~Clumping, ~Blackhole)

	// A counter so that we have some idea how it is going:
	Z := ResidueClassRing(1000);
	counter := Z!0;
	
	for i in [Minimum(Anfang, #W)..Minimum(Ende, #W)] do
		counter +:= 1;
		if counter eq Z!0 then
			printf " ... %o ... \n", i;
		end if;
		
		Enfants := ClumpingAtIndex(i, Sigma, Clumping);
		if Enfants ne { i } then
			if i in Enfants then
				printf "%o doesn't belong to sigma\n", i;
				Append(~Sigma, i);
				Append(~Clumping, Setseq(Enfants));
			elif IsEmpty(Enfants) then
				Append(~Blackhole, i);
			else
				print "shouldn't get here";
			end if;
		end if;
	end for;
end procedure;

// 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..(Rank(W)+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;

// turns a permutation to an element of the Weyl group.
permtoW := function(p)
	w := W!1;
	if #p ne Rank(W) + 1 then
		print "error";
	else
		idperm := [j : j in [1..(Rank(W) + 1)]];
		while p ne idperm do
			i := 1;
			while p[i] lt p[i+1] do
				i:= i+1;
			end while;
			x := p[i];
			p[i] := p[i+1];
			p[i+1]:= x;
			w := w * W.i;
		end while;
	end if;
	return w;
end function;

// Produces files of Sigma, f_W etc with a readable layout:
FileOutput := procedure(a,b,c)
	printf "// Description of \sigma(W) and f_W for the Coxeter group:\n";
	print W;
	
	printf "\n\n// W has order %o.\n", #W;
	
	if #c gt 0 then
		printf "\n\n// f_W is not defined on %o elements: ", #c;
		printf "u := [ \n";
		for i in [1..(#c-1)] do
			printf "%o,\n", CoxElt[c[i]];
		end for;
		printf "%o \n ];", CoxElt[c[#c]];
	end if;
	
	printf "\n\n// %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]];
	end for;
	printf "%o \n ];", CoxElt[a[#a]];
	
	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]];
		for j in b[i] do
			if j ne crit then
				printf ", %o ", CoxElt[j];
			end if;
		end for;
		printf " ],\n";
	end for;

	printf " ];";

end procedure;

// Produces files of Sigma, f_W etc with a readable layout (for type A)
FileOutputA := procedure(a,b,c)
	printf "// Description of \sigma(W) and f_W for the Coxeter group:\n";
	print W;
	
	printf "\n\n// W has order %o.\n", #W;
	
	if #c gt 0 then
		printf "\n\n// f_W is not defined on %o elements: ", #c;
		printf "u := [ \n";
		for i in [1..(#c-1)] do
			perm(CoxElt[c[i]]);
			printf ",\n";
		end for;
		perm(CoxElt[c[#c]]);
		printf " \n ];";
	end if;
	
	printf "\n\n// %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
		perm(CoxElt[a[i]]);
		printf ",\n";
	end for;
	perm(CoxElt[a[#a]]);
	printf " \n ];";
	
	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 "[ ";
		perm(CoxElt[a[i]]);
		for j in b[i] do
			if j ne crit then
				printf ", ";
				perm(CoxElt[j]);
			end if;
		end for;
		printf " ],\n";
	end for;

	printf " ];";

end procedure;

perm := procedure(w)
	p := permutation(w);
	for i in p do
		printf "%o", i;
	end for;
end procedure;