restart;with(combinat):
omega := exp((1/3)*(2*Pi*I));

P[1]:=[[1,2,3],[4,5,6],[7,8,9]]; P[2]:=[[1,4,7],[2,5,8],[3,6,9]];
P[3]:=[[1,5,9],[2,6,7],[3,4,8]]; P[4]:=[[1,6,8],[2,4,9],[3,5,7]];

rules39:={ e[1] = y[1] + y[2] + y[3] + y[4], e[2] = y[1] + omega*y[2] +
omega*y[3] + omega*y[4], e[3] = y[1] + omega^2*y[2] + omega^2*y[3] +
omega^2*y[4], e[4] = omega*y[1] + y[2] + omega^2*y[3] + omega*y[4], e[5] =
omega*y[1] + omega*y[2] + y[3] + omega^2*y[4], e[6] = omega*y[1] +
omega^2*y[2] + omega*y[3] + y[4], e[7] = omega^2*y[1] + y[2] + omega*y[3] +
omega^2*y[4], e[8] = omega^2*y[1] + omega*y[2] + omega^2*y[3] + y[4], e[9] =
omega^2*y[1] + omega^2*y[2] + y[3] + omega*y[4]
}:
rules39:=simplify(subs({y[1]=a,y[2]=b,y[3]=c,y[4]=d},rules39));

rules333:={ e[1]= z[2]+z[3]+z[4], e[2]= omega*(z[2]+z[3]+z[4]), e[3]=
omega^2*(z[2]+z[3]+z[4]), e[4]= z[2]+omega^2*z[3]+omega*z[4], e[5]= omega*z[2]
+ z[3] + omega^2*z[4], e[6]= omega^2*z[2]+ omega*z[3] + z[4], e[7]= z[2] +
omega*z[3] + omega^2*z[4], e[8]= omega*z[2] + omega^2*z[3] + z[4], e[9]=
omega^2*z[2] + z[3] + omega*z[4] };

for xx in [seq(y[i]= `+`(seq(e[P[i][1][s]],s=1..3)), i=1..4), seq(y[i]=
w^2*`+`(seq(e[P[i][2][s]],s=1..3)), i=1..4), seq(y[i]=
w*`+`(seq(e[P[i][3][s]],s=1..3)), i=1..4) ] do lprint(xx);od;

#Definition of the roots of E8 with respect to the Cartan subalgebra of sl(9)
s := 1: for ij in choose(9, 2) do beta[s] := e[ij[1]]-e[ij[2]]; beta[s+1] :=
-beta[s]; s := s+2 end do: for ijk in choose(9, 3) do beta[s] :=
e[ijk[1]]+e[ijk[2]]+e[ijk[3]]; s := s+1 end do: for ijk in choose(9, 3) do
beta[s] := -e[ijk[1]]-e[ijk[2]]-e[ijk[3]]; s := s+1 end do: s-1; s:='s':
Dad := `*`(seq(beta[i], i = 1 .. 240));
degree(Dad);
subs(rules39,Dad):e8ss:=(955004950796825236893190701774414011919935138974343129836853841/6277101735386680763835789423207666416102355444464034512896)^(-1)*factor(%);

base:=1:
for xx in seq(e8ss) do
if op(2,xx)=6 then
 base:=base*op(1,xx);
fi;
od:
B:=sort(factor(expand(base)))/4294967296;
#other fundamental invariants:
perm := proc(f,L) return
subs({x=L[1],y=L[2],z=L[3],w=L[4]},subs({a=x,b=y,c=z,d=w} ,f)); end proc;

I12 :=factor(`+`(seq(factor(subs(rules39,beta[i]^12)), i = 1 .. 240))):
I12:= I12/igcd(coeffs(I12));
z6:=1: for i from 1 to 240 do tmp:= simplify(subs(rules333, beta[i])); if
tmp<>0 then z6:=z6*tmp; fi; od: factor(z6);
z6 :=factor(`+`(seq(factor(subs(rules333,beta[i])^6), i = 1 .. 240)));
z9 :=factor(`+`(seq(factor(subs(rules333,beta[i]^9)), i = 1 .. 240)));
z12 :=factor(`+`(seq(factor(subs(rules333,beta[i]^12)), i = 1 .. 240))): z12
:= z12/igcd(coeffs(z12));

perm(I12,[a,c,b,d]);

J12:=expand(a^12+(22*b^6+(-220*c^3-220*d^3)*b^3+22*c^6-220*d^3*c^3+22*d^6)*a^6
+((-220*c^3+220*d^3)*b^6+(220*c^6-220*d^6)*b^3-220*c^6*d^3+220*c^3*d^6)*a^3+b^12
+(22*c^6+220*c^3*d^3+22*d^6)*b^6+(220*c^6*d^3+220*c^3*d^6)*b^3+c^12+22*c^6*d^6+d
^12);
perm(I12,[a,c,b,d]) - J12;

I18 :=factor(`+`(seq(simplify(subs(rules39,beta[i]^18)), i = 1 .. 240))): I18
:= (I18)/igcd(coeffs(I18));
I24 :=factor((`+`(seq(simplify(subs(rules39,beta[i]^24)), i = 1 .. 240)))):
I24 := (I24)/igcd(coeffs(I24));
I30 :=factor((`+`(seq(simplify(subs(rules39,beta[i]^30)), i = 1 .. 240)))):
I30 := (I30)/igcd(coeffs(I30));
H := CodeTools[Usage]((factor(`*`(seq(factor(subs(rules39,beta[i])), i = 1 .. 240))))):
H:=simplify((955004950796825236893190701774414011919935138974343129836853841/6277101735386680763835789423207666416102355444464034512896)^(-1) *H);
swap:=f-> subs({x=b,y=a,z=-c,w=d},subs({a=x,b=y, c=z, d=w}, f));
cycle:=f-> subs({x=b,y=-c,z=d,w=a},subs({a=x,b=y,c=z,d=w}, f));
swap(I18)-I18;
cycle(I18)-I18;
j12:=a^12+(22*b^6+(-220*c^3-220*d^3)*b^3+22*c^6-220*d^3*c^3+22*d^6)*a^6+((-220*c^3+220*d^3)*b^6+(220*c^6-220*d^6)*b^3-220*c^6*d^3+220*c^3*d^6)*a^3+b^12+(22*c^6+220*c^3*d^3+22*d^6)*b^6+(220*c^6*d^3+220*c^3*d^6)*b^3+c^12+22*c^6*d^6+d^12:
expand(j12 - perm(I12, [a,c,b,d]));

monlist:=proc(dd) local L, i12, i18, i24, i30; L := {}: for i12 from 0 to
floor(dd/12) do for i18 from 0 to floor(dd/18) do for i24 from 0 to
floor(dd/24) do for i30 from 0 to floor(dd/30) do if 12*i12 + 18*i18 + 24*i24
+ 30*i30 = dd then L := {op(L), [i12, i18, i24, i30]} end if end do end do end
do end do: L; end proc:
monlist(30);
I18;
L:=monlist(18); tmp:=`+`(seq( c[i]*I12^L[i][1]* I18^L[i][2]* I24^L[i][3]*
I30^L[i][4],i=1..nops(L) )): expand(expand(I18) - expand(tmp)):
vals:=solve({coeffs(%, [a,b,c,d])});
L:=monlist(24); tmp:=`+`(seq( c[i]*I12^L[i][1]* I18^L[i][2]* I24^L[i][3]*
I30^L[i][4],i=1..nops(L) )): expand(expand(I24) - expand(tmp)):
vals:=solve({coeffs(%, [a,b,c,d])});
L:=monlist(30); tmp:=`+`(seq( c[i]*I12^L[i][1]* I18^L[i][2]* I24^L[i][3]*
I30^L[i][4],i=1..nops(L) )): expand(expand(I30) - expand(tmp)):
vals:=solve({coeffs(%, [a,b,c,d])});
degree(H);
L:=monlist(240):nops(L); tmp:=`+`(seq( c[i]*I12^L[i][1]* I18^L[i][2]*
I24^L[i][3]* I30^L[i][4],i=1..nops(L) )): tmp2:=CodeTools[Usage](expand(H -
tmp)): vals:=solve({coeffs(tmp2, [a,b,c,d])});
roll:=rand(1000);
rndpt:=proc() {a= roll(),b= roll(), c= roll(), d= roll() }; end proc;
rndpt();
subs(rndpt(), H);

L:=monlist(240): n:=nops(L); ts:=time(): rulelist:=[seq(rndpt(), i =
1..ceil(1.1*n))]: time()-ts; ts:=time(): V:=Matrix(ceil(1.1*n),n+1): for j to
ceil(1.1*n) do printf("step %a ", j); l:=subs(rulelist[j], [I12, I18, I24,
I30]): for i to n do #printf("Case %a", i): V[j,i] :=
l[1]^L[i][1]*l[2]^L[i][2]*l[3]^L[i][3]*l[4]^L[i][4]; od; V[j,n+1]
:=subs(rulelist[j], H); od: time()-ts;
V;
CodeTools[Usage](LinearAlgebra[NullSpace](V));

Hbase:=1: for xx in seq(H) do Hbase:=Hbase*op(1,xx); od: Hbase;
factor(expand(Hbase))/2^32; disc:=%^3;
degree(Hbase);
#the gr(3,9) discriminant on semisimple elements is the cube of this base ^^

factor(`+`(seq(factor(subs(rules333,beta[i]^12)), i = 1 .. 240))):
z12 :=%/igcd(coeffs(%));

factor(`+`(seq(factor(subs(rules333,beta[i]^18)), i = 1 .. 240))):
z18 :=%/igcd(coeffs(%));

factor(`+`(seq(factor(subs(rules333,beta[i]^24)), i = 1 .. 240))):
z24 :=%/igcd(coeffs(%));

factor(`+`(seq(factor(subs(rules333,beta[i]^30)), i = 1 .. 240))):
z30 :=%/igcd(coeffs(%));

h333 := CodeTools[Usage]((factor(`*`(seq(factor(subs(rules333,beta[i])), i = 1 .. 240)))));

h333:=1: for i from 1 to 240 do
 tmp:=simplify(subs(rules333,beta[i]));
 if tmp <>0 then h333:=h333*tmp;
 fi;
od:

h333:=factor(h333)/(-1310020508637620352391208095712502073964245732475093456566329/91343852333181432387730302044767688728495783936);

degree(h333);

base:=[1,1]: for xx in seq(h333) do if op(2,xx)=6 then
base[1]:=base[1]*op(1,xx); fi; if op(2,xx)=18 then base[2]:=base[2]*op(1,xx);
fi; od: B:=sort((factor(expand(base[1]))/256)^6*
(factor(expand(base[2]))/64)^18);



#the degree 9 invariant on semisimple elements is 333 is the first factor:
z9
:=(z[3]-z[4])*(z[3]^2+z[3]*z[4]+z[4]^2)*(z[2]-z[4])*(z[2]^2+z[2]*z[4]+z[4]^2)*
(z[2]-z[3])*(z[2]^2+z[2]*z[3]+z[3]^2);
# the 333 hyperdeterminant on semisimple elements is the cube of this factor
(up to scale):

(z[2]+z[3]+z[4])*(z[2]^2-z[2]*z[3]-z[2]*z[4]+z[3]^2+2*z[3]*z[4]+z[4]^2)*(z[2]^
2-z[2]*z[3]-z[2]*z[4]+z[3]^2-z[3]*z[4]+z[4]^2)*(z[2]^2+2*z[2]*z[3]-z[2]*z[4]+z[3
]^2-z[3]*z[4]+z[4]^2)*(z[2]^2-z[2]*z[3]+2*z[2]*z[4]+z[3]^2-z[3]*z[4]+z[4]^2)*z[2
]*z[3]*z[4]; hdet333:=%^3;


katanova:=(((b^2+(2*c-d)*b+c^2-c*d+d^2)^3*(d+b+c)^3*(a^2+(-c+d)*a+(c-d)^2)^3*(
a^2+(b-d)*a+b^2+b*d+d^2)^3*(a^2+(-b-2*c)*a+b^2+b*c+c^2)^3*(b^2+(-c-d)*b+(c+d)^2)
^3*b^3*(a^2+(2*b+c)*a+b^2+b*c+c^2)^3*(a^2+(-c-2*d)*a+c^2+c*d+d^2)^3*d^3*(a^2+(-b
+c)*a+(b-c)^2)^3*(-d+a+c)^3*(-c+a+b)^3*(a^2+(-c+d)*a+c^2+c*d+d^2)^3*(a^2+(2*c+d)
*a+c^2+c*d+d^2)^3*(a^2+(-b+c)*a+b^2+b*c+c^2)^3*(a^2+(b+2*d)*a+b^2+b*d+d^2)^3*(d+
a-b)^3*(b^2+(-c+2*d)*b+c^2-c*d+d^2)^3*a^3*c^3*(a^2+(-2*b-d)*a+b^2+b*d+d^2)^3*(a^
2+(b-d)*a+(b-d)^2)^3*(b^2+(-c-d)*b+c^2-c*d+d^2)^3));
simplify(disc/katanova):numer(%);denom(%%);
perm(katanova,[a,c,b,d]);
simplify(perm(disc,[a,c,b,d])/ katanova);

#Definition of the roots of E6 with respect to the Cartan subalgebra of sl(3)sl(3)sl(3)

#the roots corresponding to the copies of sl(3)
s := 1:
for k from 0 to 2 do
beta[s]:=e[2+3*k]-e[1+3*k];beta[s+1] := -beta[s]: s := s+2:
beta[s]:=e[3+3*k]-e[1+3*k];beta[s+1] := -beta[s]: s := s+2:
beta[s]:=e[2+3*k]-e[3+3*k];beta[s+1] := -beta[s]: s := s+2:
od:
s;

#the roots corresponding to C^3C^3C^3 and its dual
for i from 1 to 3 do for j from 1 to 3 do for k from 1 to 3 do
beta[s] := e[i]+e[j+3] + e[k+6]; beta[s+1] := -beta[s]: s := s+2:
od:od:od:

dad := `*`(seq(beta[i], i = 1 .. 72));
degree(dad);

#This gives the transition between the usual cartan subalgebra (with basis x[1]..x[9] such that sum x[i]=0) and the cartan subalgebra of semi-simple 333 tensors):
omega := exp((1/3)*(2*Pi*I)):
myrules:={
e[1]= z[2]+z[3]+z[4],
e[2]= omega*(z[2]+z[3]+z[4]),
e[3]= omega^2*(z[2]+z[3]+z[4]),
e[4]= z[2]+omega^2*z[3]+omega*z[4],
e[5]= omega*z[2] + z[3] + omega^2*z[4],
e[6]= omega^2*z[2]+ omega*z[3] + z[4],
e[7]= z[2] + omega*z[3] + omega^2*z[4],
e[8]= omega*z[2] + omega^2*z[3] + z[4],
e[9]= omega^2*z[2] + z[3] + omega*z[4]};
(-2954312706550833698643/281474976710656)^(-1)*subs(myrules,dad):dd:=factor(%);
op(dd);
base:=1:
for xx in seq(dd) do
if op(2,xx)=6 then
 base:=base*op(1,xx);
fi;
od:

B:=sort(factor(expand(base)))/256;
abc:={z[2]=b, z[3]=c, z[4]= d};
subs(abc,B);