#big interpolation
with(combinat): with(LinearAlgebra): 
# increment "mylevel" 100 times - run in parallel. 
mylevel:=1:
p:=(nextprime@@mylevel)(1000);
#Definition of the roots of E7
s := 1: for ij in choose(8, 2) do beta[s] := e[ij[1]]-e[ij[2]]; beta[s+1] := -beta[s]; s := s+2 end do:
for ijkl in choose(8, 4) do beta[s] := e[ijkl[1]]+e[ijkl[2]]+e[ijkl[3]]+e[ijkl[4]]; s := s+1 end do:
s:='s':
dual_adj := product(beta[i], i = 1 .. 126):
#This gives the transition between the usual cartan subalgebra (with basis x[1]..x[8] such that sum x[i]=0) and the cartan subalgebra of semi-simple quad-vectors 
(y[1]..y[7]):
e[1] := y[1]+y[2]+y[3]+y[4]+y[5]+y[6]+y[7]:
e[2] := y[1]-y[2]+y[3]-y[4]-y[5]-y[6]+y[7]:
e[3] := y[1]+y[2]-y[3]+y[4]-y[5]-y[6]-y[7]:
e[4] := y[1]-y[2]-y[3]-y[4]+y[5]+y[6]-y[7]:
e[5] := -y[1]+y[2]+y[3]-y[4]+y[5]-y[6]-y[7]:
e[6] := -y[1]-y[2]+y[3]+y[4]-y[5]+y[6]-y[7]:
e[7] := -y[1]+y[2]-y[3]-y[4]-y[5]+y[6]+y[7]:
e[8] := -y[1]-y[2]-y[3]+y[4]+y[5]-y[6]+y[7]:

dual_adj:
factor(%):
h := %*(-1/1393796574908163946345982392040522594123776):

i:='i':
I2 := 1/288*simplify(sum(beta[i]^2, i = 1 .. 126)):
I6 := 1/6144*expand(sum(beta[i]^6, i = 1 .. 126)):
I8 :=1/16384* expand(sum(beta[i]^8, i = 1 .. 126)):
I10 := expand(sum(beta[i]^10, i = 1 .. 126)):
I10 := I10/igcd(coeffs(I10)):
I12 := expand(sum(beta[i]^12, i = 1 .. 126)):
I12 := I12/igcd(coeffs(I12)):
I14 := expand(sum(beta[i]^14, i = 1 .. 126)):
I14 := I14/igcd(coeffs(I14)):
I18 := expand(sum(beta[i]^18, i = 1 .. 126)):
I18 := I18/igcd(coeffs(I18)):
#points generator for points of I2:

roll:=rand(p):
rndpt:=proc(p)
local roll:
roll:=rand(p);
{seq(y[i] = roll(), i= 1.. 7)};
end proc:

interface(quiet =true);
 #create the partitions
L := {}:
 for i1 from 0 to 63 do
 for i2 from 0 to floor((126-2*i1)*(1/6)) do
 for i3 from 0 to floor((126-2*i1-6*i2)*(1/8)) do
 for i4 from 0 to floor((126-2*i1-6*i2-8*i3)*(1/10)) do
 for i5 from 0 to floor((126-2*i1-6*i2-8*i3-10*i4)*(1/12)) do
 for i6 from 0 to floor((126-2*i1-6*i2-8*i3-10*i4-12*i5)*(1/14)) do
 for i7 from 0 to floor((126-2*i1-6*i2-8*i3-10*i4-12*i5-14*i6)*(1/18)) do
 if 2*i1+6*i2+8*i3+10*i4+12*i5+14*i6+18*i7 = 126 then L := {op(L), [i1,i2,i3,i4,i5,i6,i7]} end if 
end do end do end do end do end do end do end do;
n:=nops(L);

ts:=time():
rulelist:=[seq(rndpt(p), i = 1..ceil(1.1*n))]:
print("made points", time()-ts);
ts:=time():
V:=Matrix(ceil(1.1*n),n+1):
for j to ceil(1.1*n) do
if j mod 1000 = 1 then printf("step %a ", j); fi;
l:=subs(rulelist[j], [I2, I6, I8, I10, I12, I14, I18]) mod p:
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]*l[5]^L[i][5]*l[6]^L[i][6]*l[7]^L[i][7] mod p; 
od;
V[j,n+1]:=- subs(rulelist[j], h ) mod p;
od:
printf("\n made big matrix in %a s\n", time()-ts);
ts:=time():
myNull:=Nullspace(V) mod p;
printf("\n Computing the null space mod %a took %a s\n", p, time()-ts);

if nops(myNull)>1 then quit; fi;
mycoeffs:=myNull[1]:
HD:= `+`(seq( mycoeffs[i]* f2^L[i][1]*f6^L[i][2]*f8^L[i][3]*f10^L[i][4]*f12^L[i][5]*f14^L[i][6]*f18^L[i][7], i= 1..n))/mycoeffs[n+1] mod p:

fname:=cat("hd48_", p);
fd:=fopen(fname, WRITE):
fprintf(fd, "COEFFS[p] := %a;\n", mycoeffs);
fclose(fd);


quit;