MAXIMAL_TERMS1 := 20000: MAXIMAL_TERMS2 := 1000: my_rank := proc (R) option remember; coxeter['rank'](R); end proc: my_cartan_matrix := proc (R) option remember; cartan_matrix(R); end proc: my_pos_roots := proc (R) option remember; pos_roots(R); end proc: ##### NUMERICAL REPRESENTATION AND REDUCTION MODULO W_P ##### reflect_weights := proc (R) option remember; local Car,i,j; Car:=my_cartan_matrix(R); [seq(omega[i]-add(Car[i,j]*omega[j],j=1..my_rank(R)),i=1..my_rank(R))]; end proc: numer_repr := proc (Car, r, w) local i, t; t:=add(omega[i],i=1..r); for i to nops(w) do t:=subs(omega[w[-i]]=Car[w[-i]],t) od; t; end proc: weyl_length := proc (Car, 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]]=Car[w[i]],t); od; res; end proc: is_minimal_repr := proc (J, Car, r, w) local i, t; t:=numer_repr(Car,r,ListTools[Reverse](w)); for i to nops(J) do if coeff(t,omega[J[i]])<0 then return false; fi; od; return true; end proc: reduce_left := proc (J, Car, r, q) local i, flag, t, res; t:=q; res:=[]; flag:=true; while flag do flag:=false; for i in J do if coeff(t,omega[i])<0 then res:=[op(res),i]; t:=subs(omega[i]=Car[i],t); flag:=true; break; fi; od; od; [res,t]; end proc: my_right_coset := proc (J, Car, r, w) local t; t:=reduce_left(J,Car,r,numer_repr (Car, r, w)); [t[1],reduce_left([$1..r], Car, r, t[2])[1]]; end proc: right_coset := proc (J, R, w) my_right_coset(J,reflect_weights(R),my_rank(R),w); end proc: my_reduce := proc (Car, r, w); my_right_coset([],Car,r,w)[2]; end proc: my_left_coset := proc (J, Car, r, w) local t; t:=my_right_coset (J, Car, r, ListTools[Reverse](w)); [my_reduce(Car,r,ListTools[Reverse](t[2])),ListTools[Reverse](t[1])]; end proc: left_coset := proc (J, R, w) my_left_coset(J,reflect_weights(R),my_rank(R),w); end proc: my_double_coset := proc (JL, JR, Car, r, w) local t; t:=my_left_coset (JR, Car, r, w);[op(my_right_coset(JL, Car, r, t[1])),t[2]]; end proc: double_coset := proc (JL, JR, R, w) my_double_coset(JL,JR,reflect_weights(R),my_rank(R),w); end proc: my_longest_elt := proc (J, R) option remember; my_right_coset(J,reflect_weights(R),my_rank(R),longest_elt(R))[1]; end proc: common_part := proc (S, R) local Car, r, T, i, j, k, flag, res; Car:=reflect_weights(R); r:=my_rank(R); T:=[seq(numer_repr(Car,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]=Car[k],T[j]); od; od; res; end proc: ##### CHOW GENERATORS ##### list_of_generators := proc (J, R, k) option remember; local Car, X, L, N, r, i, j, t, w, num; r:=my_rank(R); if k=0 then return [[[]],[add(omega[i],i=1..r)]]; fi; Car:=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,Car,r,w) then num:=numer_repr(Car,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: chow_generators := proc (J,R,k) option remember; local Car, r, L, X, i; X:=list_of_generators(J,R,k)[1]; Car:=reflect_weights(R); r:=my_rank(R); L:=[]; for i to nops(X) do L:=[op(L),my_reduce(Car,r,X[i])]; od; L; end proc: chow_dual := proc (J, R, w) local w0; w0 := my_longest_elt([$1..my_rank(R)],R); my_reduce(reflect_weights(R),my_rank(R),[op(w0),op(w),op(my_longest_elt(J,R))]); end proc: ##### PIERI FORMULA ##### pieri1 := proc (R, k, w) option remember; local Car, r, i, j, t, wop, w0, w1, res; Car:=reflect_weights(R); r:=my_rank(R); w0:=my_longest_elt([$1..r],R); wop:=my_reduce(Car,r,[op(ListTools[Reverse](w)),op(w0)]); res:=0; for i to nops(wop) do w1:=[op(wop[1..i-1]),op(wop[i+1..nops(wop)])]; if weyl_length(Car,r,w1)=nops(wop)-1 then t:=omega[k]; for j to i-1 do t:=subs(omega[wop[j]]=Car[wop[j]],t); od; res:=res+coeff(t,omega[wop[i]]) *Z[op(my_reduce(Car,r,[op(w0),op(ListTools[Reverse](w1))]))]; fi; od; res; end proc: my_indets := proc (p) local res, i, S; S:=indets(p); res:=[]; for i to nops(S) do if op(0,S[i])=Z then res:=[op(res),[seq(op(j,S[i]),j=1..nops(S[i]))]]; fi; od; res; end proc: pieri2 := proc (R, k, p) local L,i; L:=my_indets(p); if nops(L)=0 then return Z[k]*p; fi; subs(seq(Z[op(L[i])]=pieri1(R,k,L[i]),i=1..nops(L)),p) end proc: pieri := proc (R, k, n, p) local i,t; t:=p; for i to n do t:=pieri2(R,k,t); od; t; end proc: ##### DIFFERENTIAL OPERATORS ##### old_delta := proc (Car, 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)*(Car[k]-omega[k])^(m-1)*omega[k]^(i-m),m=1..i); fi; od; expand(res); end proc: critical_degree := proc (R) option remember; local r, d, n; r:=my_rank(R); n:=nops(my_pos_roots(R)); for d to n do if binomial(d+r-1,d)>MAXIMAL_TERMS1 then return d-1; fi; od; return n; end proc: my_degree := proc (r, p) local i; degree(p,{seq(omega[i],i=1..r)}); end proc: rucksack := proc (V, B, S) local A, newA, D, i, j, k, res; A:=array(0..S); newA:=array(0..S); for i from 0 to S do A[i]:=0; od; D:=array(0..nops(V),0..S); for i from 0 to S do D[0,i]:=0; od; for i to nops(V) do if V[i]=0 then for j from 0 to S do D[i,j]:=0; od; else for j from 0 to S do newA[j]:=-1; for k from 0 to min(floor(j/V[i]),B[i]) do if A[j-k*V[i]]+k*V[i]>newA[j] then newA[j]:=A[j-k*V[i]]+k*V[i]; D[i,j]:=k; fi; od; od; for j from 0 to S do A[j]:=newA[j]; od; fi; od; res:=[]; j:=S; for i from nops(V) by -1 to 1 do res:=[D[i,j],op(res)]; j:=j-V[i]*D[i,j]; od; res; end proc: my_split := proc (r, cd, p) local V, B, X, M, i, r1, r2; if op(0,p)=`^` then V:=[my_degree(r,op(1,p))]; B:=[op(2,p)]; X:=[op(1,p)]; elif op(0,p)=`*` then V:=[]; B:=[]; X:=[]; for i to nops(p) do if op(0,op(i,p))=`^` then V:=[op(V),my_degree(r,op(1,op(i,p)))]; B:=[op(B),op(2,op(i,p))]; X:=[op(X),op(1,op(i,p))]; else V:=[op(V),my_degree(r,op(i,p))]; B:=[op(B),1]; X:=[op(X),op(i,p)]; fi; od; else return [1,p]; fi; M:=rucksack(V,B,cd); r1:=1; for i to nops(V) do if M[i]>0 then r1:=r1*X[i]^M[i]; fi; od; r2:=1; for i to nops(V) do if M[i] MAXIMAL_TERMS2 then return old_delta(Car,k,expand(q)); fi; d:=my_degree(r,q); if d<=cd then return old_delta (Car,k,expand(q)); fi; if op(0,q)=`+` then return add(new_delta(Car,r,cd,k,op(i,q)),i=1..nops(q)); fi; L:=my_split(r,floor(d/2),q); L[1]:=my_expand(r,cd,L[1]); L[2]:=my_expand(r,cd,L[2]); return my_expand(r,cd,new_delta(Car,r,cd,k,L[1])*L[2])+ my_expand(r,cd,subs(omega[k]=Car[k],L[1])*new_delta(Car,r,cd,k,L[2])) end proc: c_func := proc (J, R, p) local Car, r, cd, i, j, k, L, w, w1, q, t, res; if p=0 then return 0; fi; Car:=reflect_weights(R); r:=my_rank(R); cd:=critical_degree(R); L:=chow_generators(J,R,my_degree(r,p)); w:=common_part(L,R); q:=p; for i to nops(w) do q:=new_delta(Car,r,cd,w[i],q); od; res:=0; for i to nops(L) do w1:=my_reduce(Car,r,[op(L[i]),op(w)]); t:=q; for j to nops(w1) do t:=new_delta(Car,r,cd,w1[-j],t); od; res:=res+t*Z[op(my_reduce(Car,r,L[i]))]; od; res; end proc: ##### LEVI PART COMPUTATION ##### root_sys_name := proc (type, rk); if rk=1 then return A1; fi; if substring(type,1..1)='F' then return cat(A,rk); fi; return cat(type,rk); end proc: levi_part := proc (J, Ord, R) local L, i, k, r; if Ord=[] then return []; fi; k:=0; for i to nops(Ord) do if not (Ord[i] in {op(J)}) then k:=i; break; fi; od; if k=0 then return [[Ord, R]]; fi; r:=my_rank(R); L:=substring(R,1..1); if (L='A') or (L='B') or (L='C') or (R=G2) or ((L='D') and (k>4)) or ((R=F4) and ((k=2) or (k=3))) or ((L='E') and (k>6)) then return [op(levi_part(J,Ord[1..k-1],root_sys_name(L,k-1))), op(levi_part(J,Ord[k+1..r],root_sys_name(A,r-k)))]; fi; if L='D' then if k=4 then return [op(levi_part(J,Ord[5..r],root_sys_name(A,r-4))), op(levi_part(J,[Ord[1],Ord[3],Ord[2]],A3))]; fi; if k=3 then return [op(levi_part(J,Ord[4..r],root_sys_name(A,r-3))), op(levi_part(J,[Ord[1]],A1)),op(levi_part(J,[Ord[2]],A1))]; fi; return levi_part(J,[Ord[3-k],op(Ord[3..r])],root_sys_name(A,r-1)); fi; if R=F4 then if k=1 then return levi_part(J,Ord[2..4],B3); fi; if k=4 then return levi_part(J,[Ord[3],Ord[2],Ord[1]],C3); fi; fi; if k=1 then return levi_part(J,Ord[2..r],root_sys_name(D,r-1)); fi; if k=2 then return levi_part(J,[Ord[1],op(Ord[3..r])],root_sys_name(A,r-1)); fi; if k=3 then return [op(levi_part(J,[Ord[1]],A1)), op(levi_part(J,[Ord[2],op(Ord[4..r])],root_sys_name(A,r-2)))]; fi; if k=4 then return [op(levi_part(J,[Ord[1],Ord[3]],A2)),op(levi_part(J,[Ord[2]],A1)), op(levi_part(J,Ord[5..r],root_sys_name(A,r-4)))]; fi; if k=5 then return [op(levi_part(J,[Ord[1],Ord[3],Ord[4],Ord[2]],A4)), op(levi_part(J,Ord[6..r],root_sys_name(A,r-5)))]; fi; if k=6 then return [op(levi_part(J,[Ord[2],Ord[5],Ord[4],Ord[3],Ord[1]],D5)), op(levi_part(J,Ord[7..r],root_sys_name(A,r-6)))]; fi; end proc: levi_part2 := proc (J, R) option remember; levi_part(J,[$1..my_rank(R)],R); end proc: # #### FUNDAMENTAL INVARIANTS ##### deg_fundam_invariant := proc (R, k) option remember; local L; L:=substring(R,1..1); if L='A' then return k+1; fi; if (L='B') or (L='C') then return 2*k; fi; if (L='D') then if k=2*n-1 then return c_func(J,R,subs(seq(T[i]=SUBS[i],i=1..nops(T)),p)); fi; M:=rucksack(L,B,floor(N/2)); p1:=1; d:=0; for i to nops(M) do if M[i]>0 then p1:=p1*T[i]^M[i]; fi; d:=d+L[i]*M[i]; od; p2:=1; for i to nops(M) do if M[i]=d then return c_func(J,R,subs(seq(T[i]=SUBS[i],i=1..nops(T)),p)); fi; x1:=relations(J,R,p1); G:=chow_generators(J,R,N-n); L:=generators(J,R,N-n); M:=linalg['inverse'](L[2]); V:=vector(nops(G)); for i to nops(G) do V[i]:=my_poincare_dual(J,R,x1,chow_expand_monomial(J,R,p2*L[1][i]),d); od; V:=linalg['multiply'](M,V); return add(V[i]*Z[op(chow_dual(J,R,G[i]))],i=1..nops(G)); end proc: generators := proc (J, R, n) option remember; local G, G1, L, L1, S, LP, B, M, x, i, j, k, l, rk, d, N; if n=0 then return [[1],[[1]]]; fi; G1:=chow_generators(J,R,n-1); G:=chow_generators(J,R,n); LP:=levi_part2(J,R); L:=generators(J,R,n-1); B:=[]; M:=matrix(nops(G),nops(G)); rk:=0; for i to nops(G) do for j to nops(G) do M[i,j]:=0; od; od; S:={$1..my_rank(R)} minus {op(J)}; for i to nops(S) do for j to nops(G1) do x:=pieri2(R,S[i],add(L[2][j,k]*Z[op(G1[k])],k=1..nops(G1))); for k to nops(G) do M[rk+1,k]:=coeff(x,Z[op(G[k])]) od; if linalg['rank'](M)>rk then rk:=rk+1; B:=[op(B),L[1][j]*omega[S[i]]]; fi; if rk=nops(G) then break; fi; od; if rk=nops(G) then break; fi; od; if rkrk then rk:=rk+1; B:=[op(B),P[i,j]*L1[1][l]]; if rk=nops(G) then return [B,M]; fi; fi; fi; od; fi; od; od; fi; [B,M]; end proc: add_to_mult:=proc (J,R,p) local V, L, S, SUBS, M, G, t, i, j, k; V:=my_indets(p); SUBS:=[]; for i to nops(V) do L:=generators(J,R,nops(V[i])); M:=linalg['inverse'](L[2]); G:=chow_generators(J,R,nops(V[i])); k=0; t:=my_reduce(reflect_weights(R),my_rank(R),V[i]); for j to nops(G) do if G[j]=t then k:=j; break; fi; od; if k=0 then print("WRONG GENERATOR"); return p; fi; SUBS:=[op(SUBS),add(M[k,j]*L[1][j],j=1..nops(G))]; od; subs(seq(Z[op(V[i])]=SUBS[i],i=1..nops(V)),p); end proc: c_func_inv:=proc (J,R,p) local V,LP,S,i,j; LP:=levi_part2(J,R); S:={}; for i to nops(LP) do for j to my_rank(LP[i][2]) do S:=S union {P[i,j]=fundam_invariant(LP[i][2],LP[i][1],R,j)}; od; od; subs(op(S),add_to_mult(J,R,p)); end proc: chow_expand_monomial:=proc (J,R,q) local LP, V, q1, q2, c, i, j, d; if q=0 then return 0; fi; LP:=levi_part2(J,R); V:=indets(q); c:=subs(seq(V[i]=1,i=1..nops(V)),q); for i to nops(LP) do for j to my_rank(LP[i][2]) do V:=V minus {P[i,j]}; od; od; q1:=subs(seq(V[i]=1,i=1..nops(V)),q)/c; q2:=simplify(q/q1); q2:=relations(J,R,q1)*q2; V:={$1..my_rank(R)} minus {op(J)}; for i to nops(V) do d:=degree(q2,omega[V[i]]); q2:=pieri(R,V[i],d,q2/omega[V[i]]^d); od; q2; end: chow_expand:=proc (J,R,p) local q, i, res; q:=expand(add_to_mult(J,R,p)); if op(0,q)=`+` then res:=0; for i to nops(q) do res:=res+chow_expand_monomial(J,R,op(i,q)); od; return expand(res); fi; chow_expand_monomial(J,R,q); end proc: ##### CHERN CLASSES ##### unipotent_rad := proc (J, R) option remember; local i, j, x, res, S, J2; S:=my_pos_roots(R); res:={}; J2:={$ 1..my_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: weyl_orbit := proc (J, R, v) local B, S, i, x, y, flag; B:=base(R); S:={v}; flag:=true; while flag do flag:=false; for x in S do for i in J do y:=reflect(B[i],x); if not (y in S) then flag:=true; S:=S union {y} fi; od; od; od; S; end proc: orbit_decomposition := proc (J, R, S) local T, RES, X; T:=S; RES:={}; while T<>{} do X:=weyl_orbit(J, R, T[1]); T:=T minus X; RES:=RES union {X}; od; RES; end proc: orbit_stratification := proc (J, R, S) local X, T, i; if J<>[] then X:=orbit_decomposition(J,R,S); T:={}; for i to nops(X) do T:=T union {orbit_stratification(J[1..-2],R,X[i])}; od; else T:=S; fi; T; end proc: my_flatten := proc (S) local T, X, i; if op(0,S)=`set` then T:=[]; for i to nops(S) do T:=[op(T),op(my_flatten(S[i]))]; od; else T:=[S]; fi; T; end proc: sym_funct := proc (S, R, k) local T, x, i, j; T:={}; for i to nops(S) do x:=weight_coords(S[i],R); T:=T union {add(x[j]*omega[j],j=1..my_rank(R))}; od; x:=array[0..k]; x[0]:=1; for i to k do x[i]:=0; od; for i to nops(T) do for j from k by -1 to 1 do x[j]:=x[j]+expand(x[j-1]*T[i]); od; od; x[k]; end proc: sym_functmod := proc (S, R, k, p) local T, x, i, j; T:={}; for i to nops(S) do x:=weight_coords(S[i],R); T:=T union {add(x[j]*omega[j],j=1..my_rank(R))}; od; x:=array[0..k]; x[0]:=1; for i to k do x[i]:=0; od; for i to nops(T) do for j from k by -1 to 1 do x[j]:=x[j]+expand(x[j-1]*T[i]) mod p; od; od; x[k]; end proc: chern_class := proc (J, R, k) c_func(J,R,sym_funct(unipotent_rad(J,R),R,k)); end proc: chern_classmod := proc (J, R, k, p) c_func(J,R,sym_functmod(unipotent_rad(J,R),R,k,p)) mod p; end proc: chern_class2 := proc (J, R, k) local S, T, X, i, j, m, r; r:=my_rank(R); S:=my_flatten(orbit_stratification(J, R, unipotent_rad(J, R))); T:=[]; for i to nops(S) do T:=[op(T),weight_coords(S[i],R)]; od; X:=array[0..k]; X[0]:=Z[]; for i to k do X[i]:=0; od; for i to nops(T) do for j from k by -1 to 1 do for m to r do if T[i][m]<>0 then X[j]:=X[j]+T[i][m]*pieri2(R,m,X[j-1]); fi; od; print(X); od; od; X[k]; end proc: