// 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,w> := CoxeterGroup(GrpFPCox, "A5");

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

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

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

// D5 automorphism
//a1 := hom< W-> W | [u,t,s] >;

load "A3.m";

n := Rank(W);

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

// Ideally we would like Laurent polynomials, but Magma doesn't seem to support this.
// We create a free precision ring, and this never gets us into trouble.
R<q> := LaurentSeriesRing(Integers());

// 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:
// This essentially means reversing the directions of the arrows ... 
// it is still a mystery to me why Fokko wants it around the way he does, but there is probably a good reason!
// Our W-Graph at w indexed by i is labelled by Vertices[i] , and has edges going to Edges[i], 
// which are pairs where the second paramter is the mu coefficient.
// Note: Because Fokko coded 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] := 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, arr[2]] ];
	end for;
end for;

// Elements of h are stored as sequences [  <index, poly> ] where
// index is a number between 1 and #CoxEltForms and poly is a Laurent polynomial.
// CleanHecke adds indices to make sure that every index occurs only once, and then deletes zeros.

CleanHecke := function(h)
  // first we gather like terms:
  h1 := Sort(h);
  hout := [];
  i := 1;
  while i le #h1 do
    j := i;
    while (j lt #h) and (h1[j][1] eq h1[j+1][1]) do
      j := j+1;
    end while;
    Append(~hout, <h1[i][1], &+[ h1[k][2] : k in [i..j]] >);
    i := j+1;
  end while;
  // and now we remove O's.
  for i := #hout to 1 by -1 do
    if hout[i][2] eq 0 then
      Remove(~hout, i);
    end if;
  end for;
  return hout;
end function;

// Multiplies and element in the Hecke algebra by a polynomial.
MultHecke := function(h, poly)
  return [ < h[i][1], h[i][2]*poly > : i in [1..#h] ];
end function;

// Acts on a Kazhdan-Lusztig basis
SimpleAct := function(lr, i, ind)
  if i in Vertices[ind][lr] then
    return [ < ind, q^(-1) + q >];
  else
    hout := [];
    for a in Edges[ind] do
      if i in Vertices[a[1]][lr] then
        Append(~hout, < a[1], a[2] >);
      end if;
    end for;
  end if;
  return hout;
end function;

Act := function(lr, i, h)
  hout := [];
  for a in h do
    hout cat:= MultHecke( SimpleAct(lr, i, a[1]), a[2]);
  end for;
  return CleanHecke(hout);
end function;

// This prints out an element nicely in the Cw basis.
PrintElt := procedure(h)
  h1 := CleanHecke(h);
  if h1 eq [] then
    print "0.";
  else
    printf "(%o)C(%o)", h1[1][2], CoxElt[h1[1][1]];
    for i := 2 to #h1 do
      printf " + (%o)C(%o)", h1[i][2], CoxElt[h1[i][1]];
    end for;
    printf ".";
  end if;
end procedure;

LaurentSeriesMinimum := function(f,g)
 seq1, val1 := Eltseq(f);
 seq2, val2 := Eltseq(g);

 lowerbound := Maximum(val1, val2);
 upperbound := Minimum(val1 + #seq1- 1, val2 + #seq2- 1);

 seqout := [];

 for i := lowerbound to upperbound do
   Append(~seqout, Minimum( seq1[i-val1+1], seq2[i-val2+1]) );
 end for;
 
 return q^lowerbound*(R!seqout);
end function;

HeckeMinimum := function (h1, h2);
  h1proc := CleanHecke(h1);
  h2proc := CleanHecke(h2);
  ind := [h2proc[i][1] : i in [1..#h2proc]];
  hout := [];
  for a in h1proc do
    i := Index(ind, a[1]);
    if i ne 0 then
      Append(~hout, <a[1],LaurentSeriesMinimum(R!a[2], R!h2proc[i][2])>);
    end if;
  end for;
  return CleanHecke(hout); 
end function;

CalculateBBasis := function()
  // b will be the b basis at the end:
  b := [];
  Append(~b,[<1,R!1>]);

  for i in [2..#CoxEltForms] do
    w := CoxElt[i];
    IsFirstIteration := true;
    // Remember: 1 = left, 2 = right
    for lr in [1,2] do
      for d in Vertices[i][lr] do
        if lr eq 1 then
          j := I((W.d)*w);
        else
	  j := I(w*(W.d));
        end if;
      
        if IsFirstIteration then
          btemp := Act(lr, d, b[j]);
          IsFirstIteration := false;
        else
          btemp := HeckeMinimum(btemp, Act(lr, d, b[j]));
        end if;
      end for;
    end for;
    printf ".";
    Append(~b, btemp);
  end for;
  return b;
end function;

// After some modifications to the b-basis,
// this function redoes the same algorithm to produce a (hopefully) smaller basis.
RecalculateBBasis := function(b)

  newb := [[<1,1>]];

  for i in [2..#CoxEltForms] do
    btemp := b[i];
    w := CoxElt[i];
    // Remember: 1 = left, 2 = right
    for lr in [1,2] do
      for d in Vertices[i][lr] do
        if lr eq 1 then
          j := I((W.d)*w);
        else
	  j := I(w*(W.d));
        end if;
        btemp := HeckeMinimum(btemp, Act(lr, d, newb[j]));
      end for;
    end for;
    printf ".";
    Append(~newb, btemp);
  end for;
  return newb;
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;

CriticalIndices := function(b)
  crit := [];
  for i in [1..#CoxEltForms] do
    if #b[i] gt 1 then
	Append(~crit, i);
    end if;
  end for;
  return crit;
end function;

// Calculate the T-fixed points over a point w in the Bott-Samelson resolution
Tfibre := function(BS, w)
	S := CartesianPower([0,1], #BS);
	F := [];
	for epsilon in S do
		if w eq &*[W.BS[i]^epsilon[i] : i in [1..#BS]] then
			Append(~F, epsilon);
		end if;
	end for;
	return F;
end function;

BBCellDim := function(BS, epsilon)
	dim := 0;
	prod := W!1;
	for i in [1..#BS] do
		if Length(prod) gt Length(prod* (W.BS[i])) then
			dim +:= 1;
		end if;
		prod := prod* (W.BS[i])^(epsilon[i]);
	end for;
	return dim;
end function;

BSCharacter := function(BS)
    h := [ <1, R!1> ];
    for i in BS do
		h := Act(2, i, h);
	end for;
	return h;
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;

// Produces files of Sigma, f_W etc with a readable layout:
FileOutput := procedure(b)
	printf "// Description of \sigma(W) and f_W for the Coxeter group:\n";
	print W;
	
	printf "\n\n// W has order %o.\n", #W;

	nsigmalist := [];
	for i in [1..#CoxEltForms] do
		if #b[i] gt 1 then
			Append(~nsigmalist, i);
		end if;
	end for;
	
	printf "\n\n// %o elements of W do not lie in \\sigma(W).\n// They are:\n\n", #nsigmalist;

	printf "c := [ \n";
	for i in [1..(#nsigmalist-1)] do
		printf "%o,\n", CoxElt[nsigmalist[i]];
	end for;
	printf "%o \n ];", CoxElt[nsigmalist[#nsigmalist]];
	
	printf "\n\n\n // The b-basis on the elements not in Sigma(W):\n";
	printf " // (H(w) denotes the KL-basis elt corresponding to w.) \n\n";
		
	for i in [1..#nsigmalist] do
		PrintElt(b[nsigmalist[i]]);
		printf "\n";
	end for;

end procedure;


perm := procedure(w)
	p := permutation(w);
	for i in p do
		printf "%o", i;
	end for;
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;

