numtries=1000000;(* times to try random numbers *)

(* Define lines, as vectors in I^{1,6} *)
line=Table[0,{i,27},{j,7}];

s=1;

For[i=1,i≤6,i++, line[[s,i]]=1;s++]

For[i=1,i≤6,i++, 
line[[s,7]]=2;
For[k=1,k≤6,k++, If[k≠i,line[[s,k]]=-1]];
s++
]

For[i=1,i≤6,i++, 
For[j=i+1,j≤6,j++,
line[[s,7]]=1; line[[s,i]]=-1; line[[s,j]]=-1; s++]]

(* line//MatrixForm *)

(* Define roots, as vectors in I^{1,6} *)
root=Table[0,{i,36},{j,7}];

r=1;

For[i=1,i≤6,i++, 
For[j=i+1,j≤6,j++,
root[[r,i]]=1;root[[r,j]]=-1; r++]]

For[i=1,i≤6,i++, 
For[j=i+1,j≤6,j++,
For[k=j+1,k≤6,k++,
root[[r,7]]=1;root[[r,i]]=-1;root[[r,j]]=-1; root[[r,k]]=-1;r++]]]

root[[36,7]]=2; For[i=1,i≤6,i++, root[[36,i]]=-1]

(* root//MatrixForm *)

(* Lorenzian dot product in I^{1,6} *)

mdot[x_,y_]:=x[[7]]*y[[7]]-x[[1]]*y[[1]]-x[[2]]*y[[2]]-x[[3]]*y[[3]]-x[[4]]*y[[4]]-x[[5]]*y[[5]]-x[[6]]*y[[6]];

(* paths from the center of the tree T to the endpoints.
A,B,C central branches, A0...C1 the next 6, A00...C11 the next 12 *)
ep=Table[0,{i,21},{j,12}]; (* edges-to-endpoints *)
For[i=0,i<3,i++,
ep[[7*i+1,4*i+1]]=ep[[7*i+1,4*i+2]]=ep[[7*i+1,4*i+3]]=ep[[7*i+1,4*i+4]]=1; (* A B C *)
ep[[7*i+2,4*i+1]]=ep[[7*i+2,4*i+2]]=1; (* A0 B0 C0 *)
ep[[7*i+3,4*i+3]]=ep[[7*i+3,4*i+4]]=1; (* A1 B1 C1 *)
ep[[7*i+4,4*i+1]]=1;(* A00 B00 C00 *)
ep[[7*i+5,4*i+2]]=1;(* A01 B01 C01 *)
ep[[7*i+6,4*i+3]]=1;(* A10 B10 C10 *)
ep[[7*i+7,4*i+4]]=1 (* A11 B11 C11 *)]

broot=Table[0,{i,6},{j,7}] ;(* a basis of roots of E6 *)
bnum=Table[0,{i,6}];
For[num=1,num≤numtries,num++,
For[i=1,i≤6,i++,
bnum[[i]]=RandomInteger[{1,36}];
broot[[i]]=root[[bnum[[i]]]];
];
weight=Table[0,{i,6},{j,7}];
IntForm=Table[0,{i,6},{j,6}];
For[i=1,i≤6,i++,
For[j=1,j≤6,j++,
IntForm[[i,j]]=mdot[broot[[i]],broot[[j]]]
]
];

If[Det[IntForm]==3, (* yes, these roots form a basis of E6 *)
InvForm=Inverse[IntForm];
weight=InvForm.broot;
fnum=Table[0,{i,6}];

For[i=1,i≤6,i++, fnum[[i]]=RandomInteger[{1,36}] ];
(* Print[fnum[[1]],"-",fnum[[2]],"-",fnum[[3]],"-",fnum[[4]],"-",fnum[[5]],"-",fnum[[6]]]; *)
froot=Table[0,{i,6},{j,7}];
For[i=1,i≤6,i++, froot[[i]]=root[[fnum[[i]]]] ];
l=Table[0,{e,21},{s,27},{f,6}];
For[e=1,e≤21,e++,
For[s=1,s≤27,s++,
For[f=1,f≤6,f++,
l[[e,s,f]]=ep[[e,6+f]]*mdot[froot[[f]],line[[s]]];
For[i=1,i≤6,i++,
l[[e,s,f]]=l[[e,s,f]]-ep[[e,i]]*mdot[froot[[f]],weight[[i]]]*mdot[broot[[i]],line[[s]]]; 
]
]
]
];
(* computing 21 x 27 quadratic forms, one for each of the 27 sheets *)
q=Table[0,{e,21},{s,27},{k,21}];
For[e=1,e≤21,e++,
For[s=1,s≤27,s++,
For[i=1,i<=6,i++,
For[j=1,j≤i,j++,
k=i*(i-1)/2+j;
q[[e,s,k]]=l[[e,s,i]]*l[[e,s,j]];
]
]
]
];
(* computing 21 quadratic forms *)
Q=Table[0,{e,21},{s,21}];
For[e=1,e≤21,e++,
For[s=1,s≤27,s++,
Q[[e]]+=q[[e,s]]/6
]
];
mrank=MatrixRank[Q];

If[mrank==21, (* yes, the 21 quadratic forms are linearly independent *)
Print["Found: ",Det[Q]," ",bnum[[1]],"-",bnum[[2]],"-",bnum[[3]],"-",bnum[[4]],"-",bnum[[5]],"-",bnum[[6]],"-",fnum[[1]],"-",fnum[[2]],"-",fnum[[3]],"-",fnum[[4]],"-",fnum[[5]],"-",fnum[[6]]];
q21=Table[0,{s,21},{i,6},{j,6}];
q21ranks=Table[0,{s,21}];
For[e=1,e≤21,e++,
For[i=1,i<=6,i++,
For[j=1,j≤i,j++,
k=i*(i-1)/2+j;
q21[[e,i,j]]=q21[[e,j,i]]=Q[[e,k]];
]
]
]
For[e=1,e≤21,e++,
q21ranks[[e]]=MatrixRank[q21[[e]]];
(*Print[q21[[e]]//MatrixForm];*) 
];
Print[q21ranks];
];

]
]

TimeUsed[]