needsPackage("LieAlgebraRepresentations");
debug LieAlgebraRepresentations;


-- Here we list the quantities whose factorial will be taken
-- We don't actually compute the factorials
molevFactors = (GTP) -> (
    n:=#(GTP#"weight")+1;
    p:=1;
    N1:={};
    N2:={};
    D1:={};
    D2:={};
    for k from 2 to n do (
        for i from 1 to k-1 do (
	    for j from i to k-1 do (
		N1 = append(N1,((GTP#(k,i)-i+1)-(GTP#(k-1,j)-j+1)));
		D1 = append(D1,((GTP#(k-1,i)-i+1)-(GTP#(k-1,j)-j+1)));
	    )
	);
        for i from 1 to k-1 do (
	    for j from i+1 to k do (
		N2 = append(N2,((GTP#(k,i)-i+1)-(GTP#(k,j)-j+1)-1));
		D2 = append(D2,((GTP#(k-1,i)-i+1)-(GTP#(k,j)-j+1)-1));
	    )
	);
    );
    {N1,N2,D1,D2}
)


GTPplusdelta = (GTP,i,j) -> (
    gtpA(gtpPMDeltakiEntries(GTP,1,i,j))
);


GTPminusdelta = (GTP,i,j) -> (
    gtpA(gtpPMDeltakiEntries(GTP,-1,i,j))
);


-- The minus sign comes from skipping q=i in num2
PhatoverPplus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num1 := product apply(toList(1..k-1), q -> ((GTP#(k,i)-i+1)-(GTP#(k-1,q)-q+1)+1));
    num2 := product apply(toList(1..k), q -> if q==i then 1 else ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)));
    denom1 := product apply(toList(1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)+1));
    denom2 := product apply(toList(1..k+1), q -> ((GTP#(k,i)-i+1)-(GTP#(k+1,q)-q+1)));
    (-1)*(num1*num2)/(denom1*denom2)
);


PhatoverPplusdef = (P, i, k) -> (
    Phat := GTPplusdelta(P, k,i);
    c1:=molevCoefficient(Phat);
    c2:=molevCoefficient(P);
    c1/c2
);



PhatoverPminus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num1 := product apply(toList(1..k+1), q -> ((GTP#(k,i)-i+1)-(GTP#(k+1,q)-q+1)-1));
    num2 := product apply(toList(1..k), q -> if q==i then 1 else ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)));
    denom1 := product apply(toList(1..k-1), q -> ((GTP#(k,i)-i+1)-(GTP#(k-1,q)-q+1)));
    denom2 := product apply(toList(1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)-1));
    (num1*num2)/(denom1*denom2)
);


PhatoverPminusdef = (P, i, k) -> (
    Phat := GTPminusdelta(P, k,i);
    c1:=molevCoefficient(Phat);
    c2:=molevCoefficient(P);
    c1/c2
);

end

load "Proposition_A.6_statement.m2"

-- Example 1
g = simpleLieAlgebra("A",3);
lambda = {2,0,0};
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = #Blambda;

Lplus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPPlusDeltakiEntries:=gtpPMDeltakiEntries(P,1,k,i);
      if not isValidEntryList("A",GTPPlusDeltakiEntries) then continue;
      v1=PhatoverPplus(P,i,k);
      v2=PhatoverPplusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,v1,v2,v1==v2})
)));
all(Lplus, x -> last x)

Lminus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPMinusDeltakiEntries:=gtpPMDeltakiEntries(P,-1,k,i);
      if not isValidEntryList("A",GTPMinusDeltakiEntries) then continue;
      v1=PhatoverPminus(P,i,k);
      v2=PhatoverPminusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,v1,v2,v1==v2})
)));
all(Lminus, x -> last x)




-- Example 2
g = simpleLieAlgebra("A",3);
lambda = {3,2,1};
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = #Blambda;

Lplus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPPlusDeltakiEntries:=gtpPMDeltakiEntries(P,1,k,i);
      if not isValidEntryList("A",GTPPlusDeltakiEntries) then continue;
      v1=PhatoverPplus(P,i,k);
      v2=PhatoverPplusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,v1,v2,v1==v2})
)));
all(Lplus, x -> last x)

Lminus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPMinusDeltakiEntries:=gtpPMDeltakiEntries(P,-1,k,i);
      if not isValidEntryList("A",GTPMinusDeltakiEntries) then continue;
      v1=PhatoverPminus(P,i,k);
      v2=PhatoverPminusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,v1,v2,v1==v2})
)));
all(Lminus, x -> last x)





-- Example 3
g = simpleLieAlgebra("A",5);
lambda = {1,0,0,0,1};
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = #Blambda;

Lplus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPPlusDeltakiEntries:=gtpPMDeltakiEntries(P,1,k,i);
      if not isValidEntryList("A",GTPPlusDeltakiEntries) then continue;
      v1=PhatoverPplus(P,i,k);
      v2=PhatoverPplusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,v1,v2,v1==v2})
)));
all(Lplus, x -> last x)

Lminus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPMinusDeltakiEntries:=gtpPMDeltakiEntries(P,-1,k,i);
      if not isValidEntryList("A",GTPMinusDeltakiEntries) then continue;
      v1=PhatoverPminus(P,i,k);
      v2=PhatoverPminusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,v1,v2,v1==v2})
)));
all(Lminus, x -> last x)




-- Example 4
g = simpleLieAlgebra("A",6);
lambda = {2,0,0,0,1,0};
Blambda = gtPatterns("A",dynkinToPartition("A",lambda));
N = #Blambda;

Lplus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPPlusDeltakiEntries:=gtpPMDeltakiEntries(P,1,k,i);
      if not isValidEntryList("A",GTPPlusDeltakiEntries) then continue;
      v1=PhatoverPplus(P,i,k);
      v2=PhatoverPplusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,v1,v2,v1==v2})
)));
all(Lplus, x -> last x)

Lminus = {};
for z from 0 to N-1 do (
  P = gtpA(Blambda_z);   
  for k from 1 to #lambda do (
    for i from 1 to k do (
      GTPMinusDeltakiEntries:=gtpPMDeltakiEntries(P,-1,k,i);
      if not isValidEntryList("A",GTPMinusDeltakiEntries) then continue;
      v1=PhatoverPminus(P,i,k);
      v2=PhatoverPminusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,v1,v2,v1==v2})
)));
all(Lminus, x -> last x)