# The following code is used in
# "Rationality Problem for Algebraic Tori"
# by Akinari Hoshi and Aiichi Yamasaki
# Written by Aiichi Yamasaki

DirectSumMatrixGroup:= function(l)
    local gg,gg1;
    gg:=List(l,GeneratorsOfGroup);
    if Length(Set(gg,Length))>1 then
        return fail;
    else
        gg1:=List([1..Length(gg[1])],x->DirectSumMat(List(gg,y->y[x])));
    fi;
    return Group(gg1,DirectSumMat(List(l,Identity)));
end;

DirectProductMatrixGroup:= function(l)
    local gg,gg1,o,o1,i,j,gx;
    gg:=List(l,GeneratorsOfGroup);
    gg1:=[];
    for i in [1..Length(l)] do
        o:=List(l,Identity);
        for j in gg[i] do
            o[i]:=j;
            Add(gg1,DirectSumMat(o));
        od;
    od;
    return Group(gg1,DirectSumMat(List(l,Identity)));
end;

IndmfMatrixGroup:= function(d,q,z)
    local ans;
    if d=6 and q=10 and z=1 then
        ans:=CaratMatGroupZClass(6,5517,4);
        SetName(ans,"IndmfMatrixGroup(6,10,1)");
        return ans;
    else
        return ImfMatrixGroup(d,q,z);
    fi;
end;

IndmfNumberQClasses:= function(d)
    if d=6 then
        return 10;
    else
        return ImfNumberQClasses(d);
    fi;
end;

IndmfNumberZClasses:= function(d,q)
    if d=6 and q=10 then
        return 1;
    else
        return ImfNumberZClasses(d,q);
    fi;
end;

AllImfMatrixGroups:= function(n)
    local l;
    l:=List([1..ImfNumberQClasses(n)],
      x->List([1..ImfNumberZClasses(n,x)],y->ImfMatrixGroup(n,x,y)));
    return Concatenation(l);
end;

AllIndmfMatrixGroups:= function(n)
    local l;
    l:=List([1..ImfNumberQClasses(n)],
      x->List([1..ImfNumberZClasses(n,x)],y->ImfMatrixGroup(n,x,y)));
    l:=Concatenation(l);
    if n=6 then
        l[18]:=IndmfMatrixGroup(6,10,1);
    fi;
    return l;
end;

AllDirectProductIndmfMatrixGroups:= function(l)
    local li;
    li:=List(Collected(l),
      x->UnorderedTuples(AllIndmfMatrixGroups(x[1]),x[2]));
    return List(Cartesian(li),
      x->DirectProductMatrixGroup(Concatenation(x)));
end;

SubdividedPartitions:= function(l)
    local ls;
    ls:=Cartesian(List(l,Partitions));
    return Set(ls,x->Reversed(SortedList(Concatenation(x))));
end;

PartialMatrixGroup:= function(G,l)
    local gg,gp;
    gg:=GeneratorsOfGroup(G);
    gp:=List(gg,x->x{l}{l});
    return Group(gp,IdentityMat(Length(l)));
end;

ConjugacyClassesSubgroups2:= function(g)
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    return ConjugacyClassesSubgroups(g);
end;

ConjugacyClassesSubgroupsFromPerm:= function(g)
    local iso,h,i;
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    iso:=IsomorphismPermGroup(g);
    h:=ConjugacyClassesSubgroups2(Range(iso));
    h:=List(h,Representative);
    h:=List(h,x->PreImage(iso,x));
    return h;
end;

LatticeDecompositions:= function(n)
    local d,ind,pp,p1,subgr,ld;
    ind:=[];
    for d in [1..n] do
        pp:=Partitions(d);
        p1:=List(pp,SortedList);
        if ValueOption("fromperm")=true or ValueOption("FromPerm")=true then
            subgr:=List(p1,x->
              Concatenation(List(AllDirectProductIndmfMatrixGroups(x),
              y->ConjugacyClassesSubgroupsFromPerm(y))));
        else
            subgr:=List(p1,x->
              Concatenation(List(AllDirectProductIndmfMatrixGroups(x),
              y->List(ConjugacyClassesSubgroups2(y),Representative))));
        fi;
        subgr:=List([1..NrPartitions(d)],
          x->Filtered(subgr[x],y->ForAll([1..Length(p1[x])],z->
          p1[x][z]=1 or p1[x][z]=d or
          CaratZClass(PartialMatrixGroup(y,
          [Sum([1..z-1],w->p1[x][w])+1..Sum([1..z],w->p1[x][w])]))
          in ind[p1[x][z]])));
        ld:=List(subgr,x->Set(x,CaratZClass));
        ld[NrPartitions(d)]:=Difference(ld[NrPartitions(d)],
          Union(List([1..NrPartitions(d)-1],x->ld[x])));
        ind[d]:=ld[NrPartitions(d)];
    od;
    if ValueOption("carat")=true or ValueOption("Carat")=true then
        return ld;
    else
        return List(ld,x->Set(x,Carat2CrystCat));
    fi;
end;

InverseProjection:= function(l)
    local lc,lcg,lg,ll,G,N,i,j,k,gn,gn1,gn2,o,gN,h,h1,ans,ans1;
    lc:=Collected(l);
    if ValueOption("carat")=true or ValueOption("Carat")=true then
        lcg:=List(lc,x->
          [CaratMatGroupZClass(x[1][1],x[1][2],x[1][3]),x[2]]);
    else
        lcg:=List(lc,x->
          [MatGroupZClass(x[1][1],x[1][2],x[1][3],x[1][4]),x[2]]);
    fi;
    lg:=Concatenation(List(lcg,x->List([1..x[2]],y->x[1])));
    ll:=Concatenation(List(lc,x->List([1..x[2]],y->x[1][1])));
    G:=DirectProductMatrixGroup(lg);
    gn:=[];
    for i in [1..Length(lc)] do
        if lc[i][1][1]=1 then
            gn1:=[[[-1]]];
        else
            gn1:=GeneratorsOfGroup(Normalizer(GL(lc[i][1][1],Integers),
              lcg[i][1]));
        fi;
        gn2:=[];
        for j in [1..lc[i][2]] do
            o:=List([1..lc[i][2]],x->IdentityMat(lc[i][1][1]));
            for k in gn1 do
                o[j]:=k;
                Add(gn2,DirectSumMat(o));
            od;
        od;
        gn1:=GeneratorsOfGroup(SymmetricGroup(lc[i][2]));
        gn2:=Concatenation(gn2,List(gn1,
          x->KroneckerProduct(PermutationMat(x,lc[i][2]),
          IdentityMat(lc[i][1][1]))));
        Add(gn,Group(gn2));
    od;
    N:=DirectProductMatrixGroup(gn);
    ans1:=List(ConjugacyClassesSubgroups2(G),Representative);
    ans1:=Filtered(ans1,x->
      ForAll([1..Length(lg)],y->PartialMatrixGroup(x,
      [Sum([1..y-1],z->ll[z])+1..Sum([1..y],z->ll[z])])=lg[y]));
    ans:=[];
    gN:=GeneratorsOfGroup(N);
    while ans1<>[] do
        h:=[];
        h1:=[ans1[1]];
        while h1<>[] do
            h:=Union(h,h1);
            h1:=Difference(Concatenation(List(h1,x->List(gN,y->x^y))),h);
        od;
        Add(ans,ans1[1]);
        ans1:=Difference(ans1,h);
    od;
    return ans;
end;

AllFlabbyCoflabbyZClasses:= function(n)
    local glnz,listg;
    glnz:=Concatenation(List([1..Length(cryst[n])],
      x->List([1..Length(cryst[n][x])],y->[n,x,y])));
    listg:=List(glnz,x->CaratMatGroupZClass(x[1],x[2],x[3]));
    listg:=Filtered(listg,
      x->Product(Hminus1(x))=1 and Product(H1(x))=1);
    listg:=Filtered(listg,
      x->ForAll([DerivedSubgroup(x),Centre(x),SylowSubgroup(x,2)],
      y->Product(Hminus1(y))=1 and Product(H1(y))=1));
    listg:=Filtered(listg,
      x->ForAll(List(ConjugacyClassesSubgroups2(x),Representative),
      y->Product(Hminus1(y))=1 and Product(H1(y))=1));
    if ValueOption("carat")=true or ValueOption("Carat")=true then
        return Set(listg,CaratZClass);
    else
        return Set(listg,CrystCatZClass);
    fi;
end;

AllPermutationZClasses:= function(n)
    local Sn,Snsub;
    Sn:=Group(List(GeneratorsOfGroup(SymmetricGroup(n)),
      x->PermutationMat(x,n)));
    Snsub:=List(ConjugacyClassesSubgroups2(Sn),Representative);
    if ValueOption("carat")=true or ValueOption("Carat")=true then
        return Set(Snsub,CaratZClass);
    else
        return Set(Snsub,CrystCatZClass);
    fi;
end;

MaximalGroupsID:= function(L)
    local G,O,m,m0,r,ri,ro,i,j,id,sg;
    if ValueOption("carat")=true or ValueOption("Carat")=true then
        G:=List(L,x->CaratMatGroupZClass(x[1],x[2],x[3]));
    else
        G:=List(L,x->MatGroupZClass(x[1],x[2],x[3],x[4]));
    fi;
    O:=List(G,Order);
    r:=L;
    ri:=List(r,x->Position(L,x));
    ro:=Set(ri,x->O[x]);
    m:=[];
    m0:=Filtered(ri,x->Number(ro,y->y mod O[x]=0)=1);
    while m0<>[] do
        for i in m0 do
            if ValueOption("fromperm")=true or ValueOption("FromPerm")=true then
                sg:=ConjugacyClassesSubgroupsFromPerm(G[i]);
            else
                sg:=List(ConjugacyClassesSubgroups(G[i]),Representative);
            fi;
            for j in sg do
                if Order(j) in List(ri,x->O[x]) then
                    if ValueOption("carat")=true or ValueOption("Carat")=true then
                        id:=CaratZClass(j);
                    else
                        id:=CrystCatZClass(j);
                    fi;
                    r:=Difference(r,[id]);
                fi;
            od;
        od;
        ri:=List(r,x->Position(L,x));
        ro:=Set(ri,x->O[x]);
        m:=Concatenation(m,m0);
        m0:=Filtered(ri,x->Number(ro,y->y mod O[x]=0)=1);
    od;
    return List(SortedList(m),x->L[x]);
end;