# by Viktor Petrov 2006

read("coxeter2.4v.txt"); #Stembridge
withcoxeter(); withweyl();

reflect_weights := proc (R) local C,i,j;
C:=cartan_matrix(R);
[seq(omega[i]-add(C[i,j]*omega[j],j=1..rank(R)),i=1..rank(R))];
end proc:

delta := proc (C, k, p) local res, m, c, i;
res:=0; for i to degree(p,omega[k]) do c:=coeff(p,omega[k],i); if c<>0 then res:=res+c*add(binomial(i,m)*(C[k]-omega[k])^(m-1)*omega[k]^(i-m),m=1..i); fi; od; expand(res);
end proc:

my_degree := proc (r, p) local i; degree(p,{seq(omega[i],i=1..r)}); end proc:

numer_repr := proc (C, r, w) local i, t;
t:=add(omega[i],i=1..r);
for i to nops(w) do t:=subs(omega[w[-i]]=C[w[-i]],t) od; t;
end proc:

weyl_length := proc (C, r, w) local res, i, j, t;
t:=add(omega[i],i=1..r); res := 0;
for i to nops(w) do if coeff(t,omega[w[i]])>0 then res:=res+1 else res:=res-1; fi;
t := subs(omega[w[i]]=C[w[i]],t); od; res;
end proc:

is_minimal_repr := proc (J, C, r, w) local w1, i, t;
w1:=[seq(w[-i],i=1..nops(w))]; t:=numer_repr(C,r,w1);
for i to nops(J) do if coeff(t,omega[J[i]])<0 then return false; fi; od;
return true; end proc:

common_part := proc (S, R) local C, r, T, i, j, k, flag, res;
C:=reflect_weights(R); r:=rank(R); T:=[seq(numer_repr(C,r,ListTools[Reverse](S[i])),i=1..nops(S))]; res:=[];
while true do
k:=0;
for i to r do flag:=true; for j to nops(T) do if coeff(T[j],omega[i])>0 then flag:=false; break; fi; od; if flag then k:=i; break; fi; od;
if k=0 then break; fi;
res := [op(res),k]; for j to nops(T) do T[j]:=subs(omega[k]=C[k],T[j]); od;
od;
res;
end proc:

list_of_generators := proc (J, R, k) option remember; local C, X, L, N, r, i, j, t, w, num; r:=rank(R);
if k=0 then return [[[]],[add(omega[i],i=1..r)]]; fi;
C:=reflect_weights(R);  X:=list_of_generators(J,R,k-1); N:=[]; L:=[];
for j to nops(X[1]) do for i to r do if coeff(X[2][j],omega[i])>0 then w:=[i,op(X[1][j])];
if is_minimal_repr(J,C,r,w) then num:=numer_repr(C,r,w);
if not (num in {op(N)}) then L:=[op(L),w]; N:=[op(N),num]; fi; fi; fi; od; od; [L,N];
end proc:

c_func := proc (J, R, p) local C, r, cd, i, j, k, L, w, w1, q, t, res;
if p=0 then return 0; fi;
C:=reflect_weights(R); L:=list_of_generators(J,R,my_degree(rank(R),p))[1];
w:=common_part(L,R); q:=p; for i to nops(w) do q:=delta(C,w[i],q); od; res:=0;
for i to nops(L) do w1:=reduce([op(L[i]),op(w)],R); t:=q; for j to nops(w1) do t:=delta(C,w1[-j],t); od; res:=res+t*Z[op(reduce(L[i],R))]; od; res;
end proc:

unipotent_rad := proc (J, R) local i, j, x, res, S, J2;
S:=pos_roots(R); res:={}; J2:={$ 1..rank(R)} minus {op(J)};
for i to nops(S) do x:=root_coords (S[i], R); for j to nops(J2) do if x[J2[j]]>0 then res:=res union {S[i]}; break; fi; od; od; res;
end proc:

sym_funct := proc (S, k) local x, i, j;
x:=array[0..k];
x[0]:=1; for i to k do x[i]:=0; od;
for i to nops(S) do for j from k by -1 to 1 do x[j]:=x[j]+expand(x[j-1]*S[i]); od;
od;
x[k];
end proc:

chern_class := proc (J, R, k) local S, T, x, i, j;
S:=unipotent_rad(J,R); T:={};
for i to nops(S) do x:=weight_coords(S[i],R); T:=T union {add(x[j]*omega[j],j=1..rank(R))}; od;
c_func(J,R,sym_funct(T,k));
end proc:

# Example which computes the 10-th Chern class of the tangent bundle of E8/P8
#chern_class([1,2,3,4,5,6,7],E8,10);