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))
);


N1hatoverN1plus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num:=1;
    if i<k then (
	num = product apply(toList(i..k-1), q -> ((GTP#(k,i)-i+1)-(GTP#(k-1,q)-q+1)+1))
    );
    denom := product apply(toList(1..i), p -> ((GTP#(k+1,p)-p+1)-(GTP#(k,i)-i+1)));
    num/denom
);


N1hatoverN1plusdef = (P, i, k) -> (
    Phat := GTPplusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_0, t -> t!))/(product apply(M2_0, t -> t!))
);


D1hatoverD1plus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num:=1;
    if i<k then (
	num = product apply(toList(i+1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)+1))
    );
    denom:=1;
    if i>1 then (
        denom = product apply(toList(1..i-1), p -> ((GTP#(k,p)-p+1)-(GTP#(k,i)-i+1)))
    );
    num/denom
);


D1hatoverD1plusdef = (P, i, k) -> (
    Phat := GTPplusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_2, t -> t!))/(product apply(M2_2, t -> t!))
);


N2hatoverN2plus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num:=1;
    if i<k then (
	num = product apply(toList(i+1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)))
    );
    denom:=1;
    if i>1 then (
        denom = product apply(toList(1..i-1), p -> ((GTP#(k,p)-p+1)-(GTP#(k,i)-i+1)-1))
    );
    num/denom
);


N2hatoverN2plusdef = (P, i, k) -> (
    Phat := GTPplusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_1, t -> t!))/(product apply(M2_1, t -> t!))
);


D2hatoverD2plus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num := product apply(toList(i+1..k+1), q -> ((GTP#(k,i)-i+1)-(GTP#(k+1,q)-q+1)));
    denom:=1;
    if i>1 then (
        denom = product apply(toList(1..i-1), p -> ((GTP#(k-1,p)-p+1)-(GTP#(k,i)-i+1)-1))
    );
    num/denom
);


D2hatoverD2plusdef = (P, i, k) -> (
    Phat := GTPplusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_3, t -> t!))/(product apply(M2_3, t -> t!))
);



N1hatoverN1minus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num := product apply(toList(1..i), p -> ((GTP#(k+1,p)-p+1)-(GTP#(k,i)-i+1)+1));
    denom := product apply(toList(i..k-1), q -> ((GTP#(k,i)-i+1)-(GTP#(k-1,q)-q+1)));
    num/denom
);


N1hatoverN1minusdef = (P, i, k) -> (
    Phat := GTPminusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_0, t -> t!))/(product apply(M2_0, t -> t!))
);


D1hatoverD1minus= (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num:=product apply(toList(1..i-1), p -> ((GTP#(k,p)-p+1)-(GTP#(k,i)-i+1))+1);
    denom:=product apply(toList(i+1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)));
    num/denom
);


D1hatoverD1minusdef = (P, i, k) -> (
    Phat := GTPminusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_2, t -> t!))/(product apply(M2_2, t -> t!))
);


N2hatoverN2minus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num:=product apply(toList(1..i-1), p -> ((GTP#(k,p)-p+1)-(GTP#(k,i)-i+1)));
    denom:=product apply(toList(i+1..k), q -> ((GTP#(k,i)-i+1)-(GTP#(k,q)-q+1)-1));
    num/denom
);


N2hatoverN2minusdef = (P, i, k) -> (
    Phat := GTPminusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_1, t -> t!))/(product apply(M2_1, t -> t!))
);

D2hatoverD2minus = (GTP, i,k) -> (
    n:=#(GTP#"weight")+1;
    num:=product apply(toList(1..i-1), p -> ((GTP#(k-1,p)-p+1)-(GTP#(k,i)-i+1)));
    denom := product apply(toList(i+1..k+1), q -> ((GTP#(k,i)-i+1)-(GTP#(k+1,q)-q+1)-1));
    num/denom
);


D2hatoverD2minusdef = (P, i, k) -> (
    Phat := GTPminusdelta(P, k,i);
    M1:=molevFactors(Phat);
    M2:=molevFactors(P);
    (product apply(M1_3, t -> t!))/(product apply(M2_3, t -> t!))
);

end

load "Proposition_A.6_proof.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;
      -- Check the formulas against the definitions
      b1=N1hatoverN1plus(P,i,k)==N1hatoverN1plusdef(P,i,k);
      b2=D1hatoverD1plus(P,i,k)==D1hatoverD1plusdef(P,i,k);
      b3=N2hatoverN2plus(P,i,k)==N2hatoverN2plusdef(P,i,k);
      b4=D2hatoverD2plus(P,i,k)==D2hatoverD2plusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions      
      b1=N1hatoverN1minus(P,i,k)==N1hatoverN1minusdef(P,i,k);
      b2=D1hatoverD1minus(P,i,k)==D1hatoverD1minusdef(P,i,k);
      b3=N2hatoverN2minus(P,i,k)==N2hatoverN2minusdef(P,i,k);
      b4=D2hatoverD2minus(P,i,k)==D2hatoverD2minusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions
      b1=N1hatoverN1plus(P,i,k)==N1hatoverN1plusdef(P,i,k);
      b2=D1hatoverD1plus(P,i,k)==D1hatoverD1plusdef(P,i,k);
      b3=N2hatoverN2plus(P,i,k)==N2hatoverN2plusdef(P,i,k);
      b4=D2hatoverD2plus(P,i,k)==D2hatoverD2plusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions      
      b1=N1hatoverN1minus(P,i,k)==N1hatoverN1minusdef(P,i,k);
      b2=D1hatoverD1minus(P,i,k)==D1hatoverD1minusdef(P,i,k);
      b3=N2hatoverN2minus(P,i,k)==N2hatoverN2minusdef(P,i,k);
      b4=D2hatoverD2minus(P,i,k)==D2hatoverD2minusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions
      b1=N1hatoverN1plus(P,i,k)==N1hatoverN1plusdef(P,i,k);
      b2=D1hatoverD1plus(P,i,k)==D1hatoverD1plusdef(P,i,k);
      b3=N2hatoverN2plus(P,i,k)==N2hatoverN2plusdef(P,i,k);
      b4=D2hatoverD2plus(P,i,k)==D2hatoverD2plusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions      
      b1=N1hatoverN1minus(P,i,k)==N1hatoverN1minusdef(P,i,k);
      b2=D1hatoverD1minus(P,i,k)==D1hatoverD1minusdef(P,i,k);
      b3=N2hatoverN2minus(P,i,k)==N2hatoverN2minusdef(P,i,k);
      b4=D2hatoverD2minus(P,i,k)==D2hatoverD2minusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions
      b1=N1hatoverN1plus(P,i,k)==N1hatoverN1plusdef(P,i,k);
      b2=D1hatoverD1plus(P,i,k)==D1hatoverD1plusdef(P,i,k);
      b3=N2hatoverN2plus(P,i,k)==N2hatoverN2plusdef(P,i,k);
      b4=D2hatoverD2plus(P,i,k)==D2hatoverD2plusdef(P,i,k);
      Lplus = append(Lplus,{z,k,i,all({b1,b2,b3,b4})})
)));
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;
      -- Check the formulas against the definitions      
      b1=N1hatoverN1minus(P,i,k)==N1hatoverN1minusdef(P,i,k);
      b2=D1hatoverD1minus(P,i,k)==D1hatoverD1minusdef(P,i,k);
      b3=N2hatoverN2minus(P,i,k)==N2hatoverN2minusdef(P,i,k);
      b4=D2hatoverD2minus(P,i,k)==D2hatoverD2minusdef(P,i,k);
      Lminus = append(Lminus,{z,k,i,all({b1,b2,b3,b4})})
)));
all(Lminus, x -> last x)