--first choose your n

n=10

--load into macaulay a file containing hypertree divisors Di, as well as a collection of n! row vectors encoding the symmeric group. (we generated the symmetric group in matlab and imported it to macaulay.)

--Define SG to be your set of vectors encoding the symmeric group.

--create a matrix with te Di's as entries. Note that the number of entries in MM must be altered for each n. 

--The point of this step is to write a function D with D(i)=Di (ith hypertree divisor) without having to run the entire hypertree code. We can then easily loop over all divisor classes.


MM=matrix{{D0,D1,D2,D3,D4,D5,D6,D7,D8,D9,D10,D11,D12,D13,D14,D15,D16,D17,D18,D19,D20,D21,D22,D23,D24,D25,D26,D27,D28,D29,D30,D31,D32,D33,D34,D35,D36,D37,D38,D39,D40,D41,D42,D43,D44,D45,D46,D47,D48,D49,D50,D51,D52,D53,D54,D55,D56,D57,D58,D59,D60,D61,D62,D63,D64,D65,D66,D67,D68,D69,D70,D71,D72,D73,D74,D75,D76,D77,D78,D79,D80,D81,D82,D83,D84,D85,D86,D87,D88,D89,D90,D91,D92,D93,D94,D95}};

D = i-> MM_(0,i);


--computing symmetries amounts to substitutions based on each row vector in SG.


gg = (i,v)-> if i< n then x_(v_(0,i)) else 0;

sb = (P,k)-> sub(P,{x_0=>gg(0,k), x_1=>gg(1,k), x_2=>gg(2,k), x_3=>gg(3,k), x_4=>gg(4,k), x_5=>gg(5,k), x_6=>gg(6,k),x_7=>gg(7,k),x_8=>gg(8,k)});



--(Problem 1) determine size of symmetry group of divisor. for divisor equation p and set SG of all possible permutations, sz(p) returns permutations that preserve p.


ss = (p,l) -> if p==sb(p,l) then true else false;
sy = (p,L) -> select( toList( set{p}**set L), ss);

aa = (x,y) ->y;
sz = p-> apply( sy(p,SG) , aa);

scan(0..nn(n), i-> <<"D"<<i<<" has symmetry group order "<<num(sz(D(i))) )



--(Problem 2) determine when pairs of divisors give the same class for some permutations of {1,..n}

--note that this process can be curtailed by first inspecting symmetry groups sizes. Evidently Di and Dj are permutations of eachother only if they have the same symmetry group.

jj = p -> set(  apply(toList( set {p} ** set( SG ) ), sb) );

ee = (p,q) -> if toList (jj(p)*jj(q)) == {} then false else true;

eee=(i,j)->(i,j,ee(D(i),D(j)));

ttt=(i,j)->if i<j then true else false;
K= select( toList ( (set for i from 0 to nn(n) list i)** (set for i from 0 to nn(n) list i) ) , ttt);


--next line will output (i,j,true/false); (i,j,true) indicates ht divisors Di and Dj are permutations of eachother; false indicates they are not.

apply(K,eee)