restart; with(LinearAlgebra); interface(rtablesize = 40);
#Katanova_Invariants (computed by another method and copy-pasted)

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:

j18:=a^18+(-17*b^6+(170*c^3+170*d^3)*b^3-17*c^6+170*d^3*c^3-17*d^6)*a^12+((1870*c^3-1870*d^3)*b^6+(-1870*c^6+1870*d^6)*b^3+1870*c^6*d^3-1870*c^3*d^6)*a^9+(-17*b^12+(1870*c^3+1870*d^3)*b^9+(-7854*c^6-7854*d^6)*b^6+(1870*c^9+1870*d^9)*b^3-17*c^12+1870*c^9*d^3-7854*c^6*d^6+1870*c^3*d^9-17*d^12)*a^6+((170*c^3-170*d^3)*b^12+(-1870*c^6+1870*d^6)*b^9+(1870*c^9-1870*d^9)*b^6+(-170*c^12+170*d^12)*b^3+170*c^12*d^3-1870*c^9*d^6+1870*c^6*d^9-170*c^3*d^12)*a^3+b^18+(-17*c^6-170*c^3*d^3-17*d^6)*b^12+(-1870*c^6*d^3-1870*c^3*d^6)*b^9+(-17*c^12-1870*c^9*d^3-7854*c^6*d^6-1870*c^3*d^9-17*d^12)*b^6+(-170*c^12*d^3-1870*c^9*d^6-1870*c^6*d^9-170*c^3*d^12)*b^3+c^18-17*c^12*d^6-17*c^6*d^12+d^18:

j24:=111*a^24+(506*b^6+(-5060*c^3-5060*d^3)*b^3+506*c^6-5060*d^3*c^3+506*d^6)*a^18+((-206448*c^3+206448*d^3)*b^6+(206448*c^6-206448*d^6)*b^3-206448*c^6*d^3+206448*c^3*d^6)*a^15+(10166*b^12+(-1118260*c^3-1118260*d^3)*b^9+(4696692*c^6+4696692*d^6)*b^6+(-1118260*c^9-1118260*d^9)*b^3+10166*c^12-1118260*c^9*d^3+4696692*c^6*d^6-1118260*c^3*d^9+10166*d^12)*a^12-1118260*(b-d)*(b^6+(-10*c^3-10*d^3)*b^3+c^6-10*d^3*c^3+d^6)*(b-c)*(c-d)*(b^2+b*c+c^2)*(c^2+c*d+d^2)*(b^2+b*d+d^2)*a^9+(506*b^18+(-206448*c^3-206448*d^3)*b^15+(4696692*c^6+4696692*d^6)*b^12+(-12300860*c^9-12300860*d^9)*b^9+(4696692*c^12+4696692*d^12)*b^6+(-206448*c^15-206448*d^15)*b^3+506*c^18-206448*c^15*d^3+4696692*c^12*d^6-12300860*c^9*d^9+4696692*c^6*d^12-206448*c^3*d^15+506*d^18)*a^6+((-5060*c^3+5060*d^3)*b^18+(206448*c^6-206448*d^6)*b^15+(-1118260*c^9+1118260*d^9)*b^12+(1118260*c^12-1118260*d^12)*b^9+(-206448*c^15+206448*d^15)*b^6+(5060*c^18-5060*d^18)*b^3-5060*c^18*d^3+206448*c^15*d^6-1118260*c^12*d^9+1118260*c^9*d^12-206448*c^6*d^15+5060*c^3*d^18)*a^3+111*b^24+(506*c^6+5060*c^3*d^3+506*d^6)*b^18+(206448*c^6*d^3+206448*c^3*d^6)*b^15+(10166*c^12+1118260*c^9*d^3+4696692*c^6*d^6+1118260*c^3*d^9+10166*d^12)*b^12+(1118260*c^12*d^3+12300860*c^9*d^6+12300860*c^6*d^9+1118260*c^3*d^12)*b^9+(506*c^18+206448*c^15*d^3+4696692*c^12*d^6+12300860*c^9*d^9+4696692*c^6*d^12+206448*c^3*d^15+506*d^18)*b^6+(5060*c^18*d^3+206448*c^15*d^6+1118260*c^12*d^9+1118260*c^9*d^12+206448*c^6*d^15+5060*c^3*d^18)*b^3+111*c^24+506*c^18*d^6+10166*c^12*d^12+506*c^6*d^18+111*d^24:

j30:=584*a^30+(-435*b^6+(4350*c^3+4350*d^3)*b^3-435*c^6+4350*d^3*c^3-435*d^6)*a^24+((440220*c^3-440220*d^3)*b^6+(-440220*c^6+440220*d^6)*b^3+440220*c^6*d^3-440220*c^3*d^6)*a^21+(-63365*b^12+(6970150*c^3+6970150*d^3)*b^9+(-29274630*c^6-29274630*d^6)*b^6+(6970150*c^9+6970150*d^9)*b^3-63365*c^12+6970150*c^9*d^3-29274630*c^6*d^6+6970150*c^3*d^9-63365*d^12)*a^18+25852920*(b-d)*(b^6+(-10*c^3-10*d^3)*b^3+c^6-10*d^3*c^3+d^6)*(b-c)*(c-d)*(b^2+b*c+c^2)*(c^2+c*d+d^2)*(b^2+b*d+d^2)*a^15+(-63365*b^18+(25852920*c^3+25852920*d^3)*b^15+(-588153930*c^6-588153930*d^6)*b^12+(1540403150*c^9+1540403150*d^9)*b^9+(-588153930*c^12-588153930*d^12)*b^6+(25852920*c^15+25852920*d^15)*b^3-63365*c^18+25852920*c^15*d^3-588153930*c^12*d^6+1540403150*c^9*d^9-588153930*c^6*d^12+25852920*c^3*d^15-63365*d^18)*a^12+6970150*(b^12+(-199/5*c^3-199/5*d^3)*b^9+(906/5*c^6+702/5*d^3*c^3+906/5*d^6)*b^6+(-199/5*c^9+702/5*c^6*d^3+702/5*c^3*d^6-199/5*d^9)*b^3+c^12-199/5*c^9*d^3+906/5*c^6*d^6-199/5*c^3*d^9+d^12)*(b-d)*(b-c)*(c-d)*(b^2+b*c+c^2)*(c^2+c*d+d^2)*(b^2+b*d+d^2)*a^9+(-435*b^24+(440220*c^3+440220*d^3)*b^21+(-29274630*c^6-29274630*d^6)*b^18+(284382120*c^9+284382120*d^9)*b^15+(-588153930*c^12-588153930*d^12)*b^12+(284382120*c^15+284382120*d^15)*b^9+(-29274630*c^18-29274630*d^18)*b^6+(440220*c^21+440220*d^21)*b^3-435*c^24+440220*c^21*d^3-29274630*c^18*d^6+284382120*c^15*d^9-588153930*c^12*d^12+284382120*c^9*d^15-29274630*c^6*d^18+440220*c^3*d^21-435*d^24)*a^6+((4350*c^3-4350*d^3)*b^24+(-440220*c^6+440220*d^6)*b^21+(6970150*c^9-6970150*d^9)*b^18+(-25852920*c^12+25852920*d^12)*b^15+(25852920*c^15-25852920*d^15)*b^12+(-6970150*c^18+6970150*d^18)*b^9+(440220*c^21-440220*d^21)*b^6+(-4350*c^24+4350*d^24)*b^3+4350*d^3*c^24-440220*d^6*c^21+6970150*d^9*c^18-25852920*d^12*c^15+25852920*d^15*c^12-6970150*d^18*c^9+440220*d^21*c^6-4350*d^24*c^3)*a^3+584*b^30+(-435*c^6-4350*c^3*d^3-435*d^6)*b^24+(-440220*c^6*d^3-440220*c^3*d^6)*b^21+(-63365*c^12-6970150*c^9*d^3-29274630*c^6*d^6-6970150*c^3*d^9-63365*d^12)*b^18+(-25852920*c^12*d^3-284382120*c^9*d^6-284382120*c^6*d^9-25852920*c^3*d^12)*b^15+(-63365*c^18-25852920*c^15*d^3-588153930*c^12*d^6-1540403150*c^9*d^9-588153930*c^6*d^12-25852920*c^3*d^15-63365*d^18)*b^12+(-6970150*c^18*d^3-284382120*c^15*d^6-1540403150*c^12*d^9-1540403150*c^9*d^12-284382120*c^6*d^15-6970150*c^3*d^18)*b^9+(-435*c^24-440220*c^21*d^3-29274630*c^18*d^6-284382120*c^15*d^9-588153930*c^12*d^12-284382120*c^9*d^15-29274630*c^6*d^18-440220*c^3*d^21-435*d^24)*b^6+(-4350*c^24*d^3-440220*c^21*d^6-6970150*c^18*d^9-25852920*c^15*d^12-25852920*c^12*d^15-6970150*c^9*d^18-440220*c^6*d^21-4350*c^3*d^24)*b^3+584*c^30-435*c^24*d^6-63365*c^18*d^12-63365*c^12*d^18-435*c^6*d^24+584*d^30:

#check invariance
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));
for xx in {j12,j18,j24,j30} do expand(swap(xx)-xx) od;
for xx in {j12,j18,j24,j30} do expand(cycle(xx)-xx) od;


#Create partition

L := []; for i from 0 to 10 do for j from 0 to 6 do for k from 0 to 5 do for l from 0 to 4 do if 12*i+18*j+24*k+30*l = 120 then L := [op(L), [i, j, k, l]] end if end do end do end do end do;

L;
#Generate a vector of monomials
for i to 28 do V[i] := j12^L[i][1]*j18^L[i][2]*j24^L[i][3]*j30^L[i][4] end do;

for i to 28 do F[i] := f[12]^L[i][1]*f[18]^L[i][2]*f[24]^L[i][3]*f[30]^L[i][4] end do;

#Create the matrix of 28 monomials evaluated on 4*pp points of the dual
roll:=rand(-10..10);
pp:=9:
M := Matrix(4*pp, 28):
for i to pp do
p2:= roll();
p3:= roll();
p4:= roll();
for j to 28 do
	M[i, j] := expand(subs([a = 0, b = p2, c = p3, d = p4], V[j])) end do:
	p2:= roll(); p3:= roll(); p4:= roll();
	for j to 28 do M[i+pp,j]:= subs([b = 0, a = p2, c = p3, d = p4], V[j]) end do:
	p2:= roll(); p3:= roll(); p4:= roll();
	for j to 28 do M[i+2*pp,j] := subs([c = 0, a = p2, b =p3, d = p4], V[j]) end do:
	p2:= roll(); p3:= roll(); p4:= roll();
	for j to 28 do M[i+3*pp,j] := subs([d = 0, a = p2, b = p3, c = p4], V[j]) end do:
od:
Rank(M);

#this vector of rational numbers gives the relation among the fundamental invariants
roll20 :=rand(10^20);
nspace:=seq(NullSpace(M)); hf:=`+`(seq(nspace[i]*F[i], i = 1..28))*roll20();
inspace := ilcm(seq(denom(nspace[i]), i=1..28))*nspace;

HF:=`+`(seq(inspace[i]*F[i], i = 1..28));

#practice what to do mod several primes.
`mod`:=mods;

N:=30: P:=[seq((nextprime@@i)(300), i = 1..N) ];
C:=Matrix([seq([seq(nspace[j] mod P[i], j=1..28)], i = 1..N)]);
`*`(seq(P));
n:=28;
chf:=[seq(chrem([seq(C[i,j], i=1..N-1)], P[1..N-1]), j=1..n)];
M:=`*`(seq(P[1..N-1])); nv :=[seq(iratrecon(chf[i], M),i = 1..n)];

mypoly:=`+`(seq(nv[i]*F[i], i = 1..28));
Matrix(n,1,nv) - Matrix(n,1,nspace);
hf:=`+`(seq(nspace[i]*F[i], i = 1..28)): sort(%);
mypoly - hf;
quit;