load "FermatTwists.m";
load "KleinTwists.m";

n := 4;

Fsimplices:=AssociativeArray();
for i in [1..#Ftwists] do
    G0,H,H0 := groups(Ftwists[i],Fphi);
    for j in [1,2] do CC:=CorrelationSimplex(G0,H,H0,j,n,Fexp); keys:=Sort([v:v in Keys(CC)]); Fsimplices[<i,j>]:=[<v,CC[v]>:v in keys]; end for;
    print i; 
end for;

Ksimplices:=AssociativeArray();
for i in [1..#Ktwists] do
    G0,H,H0 := groups (Ktwists[i],Kphi);
    for j in [1,2] do CC:=CorrelationSimplex(G0,H,H0,j,n,Kexp); keys:=Sort([v:v in Keys(CC)]); Ksimplices[<i,j>]:=[<v,CC[v]>:v in keys]; end for;
    print i; 
end for;

Fdists := {Fsimplices[v]:v in Keys(Fsimplices)};
print #Fdists, "distinct Fermat twist distributions with n =", n;

Kdists := {Ksimplices[v]:v in Keys(Ksimplices)};
print #Kdists, "distinct Klein twist distributions with n =", n;

print #(Fdists join Kdists), "distinct twists distributions in total with n =", n;

print #(Kdists meet Fdists), "distinct Klein twist distributions are also Fermat twist distributions with n =", n;

STdict, STgroups := STdictionary (<<Ftwists,Fphi,Fexp>,<Ktwists,Kphi,Kexp>>);
printf "STgroups contains %o distinct distributions indexed by Stdict with key <M101,M030,M202>\n", #STgroups;

V := [<1,0,1>,<0,3,0>,<2,0,2>,<2,0,0>,<4,0,0>,<0,1,0>,<0,2,0>,<0,0,2>,<0,0,4>];
iV := [<0,3,0>,<2,0,0>,<4,0,0>,<0,1,0>,<0,2,0>,<0,0,2>,<0,0,4>]; 

function FKstats(prefix,twist,j)
    M:=AssociativeArray();
    S := Split(Read(Sprintf("%o%o_%o_moments.txt",prefix,twist[1],j)),"\n");
    r := Split(S[#S]," ");
    z1 := Round(48*eval(r[10]))/48;
    z2 := [Round(48*a)/48: a in eval(r[11])];
    z3 := Round(48*eval(r[12]))/48;
    return [eval(r[i]):i in [1..9]],z1,z2,z3;
end function;

function FKexps(twist,j,phi,exp)
    G,H,H0 := groups(twist,phi);
    M:=AssociativeArray();
    for v in V do M[v] := vExp (G,H,H0,j,v,exp); end for;
    z1 := z1Exp(H,H0,j);
    z2 := z2Exp(H,H0,j);
    z3 := z3Exp(H,H0,j);
    return [vExp (G,H,H0,j,v,exp):v in V],z1,z2,z3;
end function;

function FKverify_stats(twists,phi,exp);
    assert #twists in [14,59];
    prefix := #twists eq 59 select "../Moments/FermatTwist" else "../Moments/KleinTwist";
    for r in twists do
        for j in [1,2] do
            if j eq 1 then
                G,H,H0 := groups(r,phi);
                if H eq H0 then continue; end if;
            end if;
            Ms,z1s,z2s,z3s := FKstats(prefix,r,j);
            Me,z1e,z2e,z3e := FKexps(r,j,phi,exp);
            assert z1s eq z1e and z2s eq z2e and z3s eq z3e;
            for i:=1 to #V do
                if Abs(1.0-Ms[i]/Me[i]) ge 0.01 then print r[1],j; end if;
            end for;
        end for;
    end for;
    return true;
end function;

STvecs :=[<vExp(a[1],a[2],a[3],a[4],v,a[5]):v in V>:a in STgroups];
STivecs :=[<vExp(a[1],a[2],a[3],a[4],v,a[5]):v in iV>:a in STgroups];

print "Pairs of ST distributions with the same independent coefficient distributions:";
for i in [1..#STgroups] do for j in [i+1..#STgroups] do
    if STivecs[i] eq STivecs[j] then print i,j; end if;
end for; end for;


print "Pairs of ST distributions that differ by less than one percent:";
for i in [1..#STgroups] do for j in [i+1..#STgroups] do
    if #[k:k in [1..#V]|Abs(1-STvecs[i][k]/STvecs[j][k]) gt 0.01] eq 0 then print i,j; end if;
end for; end for;
