# The following code is used in
# "Multiplicative Invariant Fields of Dimension <=6"
# by Akinari Hoshi, Ming-Chang Kang and Aiichi Yamasaki
# Written by Aiichi Yamasaki

NormalSubgroups2:= function(g)
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    return NormalSubgroups(g);
end;

DoubleCosets2:= function(g,u,v)
    Reset(GlobalMersenneTwister);
    Reset(GlobalRandomSource);
    return DoubleCosets(g,u,v);
end;

OuterConjugateSubgroups:= function(N,H)
    local n,hc,hc1,hc2;
    n:=GeneratorsOfGroup(N);
    hc:=[];
    hc1:=[H];
    while hc1<>[] do
        Append(hc,hc1);
        hc2:=Flat(List(hc1,x->List(n,y->x^y)));
        hc1:=Difference(hc2,hc);
    od;
    return hc;
end;

AllSubdirectProducts:= function(G1,G2)
    local G1g,Ng1,Ng2,N1,N2,nn1,nn2,Nnn1,Nnn2,n,nn,nc,ans,n1,n2,epi1,epi2,iso,aut,cc1,cc2,c1,c2,d,dr,f;
    G1g:=GeneratorsOfGroup(G1);
    if Rank(Identity(G1))=1 then
        Ng1:=[[[-1]]];
    else
        Ng1:=Unique(GeneratorsOfGroup(NormalizerInGLnZ(G1)));
    fi;
    if Rank(Identity(G2))=1 then
        Ng2:=[[[-1]]];
    else
        Ng2:=Unique(GeneratorsOfGroup(NormalizerInGLnZ(G2)));
    fi;
    N1:=NormalSubgroups2(G1);
    nn1:=[];  Nnn1:=[];
    while N1<>[] do
        n:=N1[1];
        nn:=Unique(GeneratorsOfGroup(Normalizer(Group(Ng1),n)));
        Add(nn1,n);
        Add(Nnn1,nn);
        nc:=OuterConjugateSubgroups(Group(Ng1),n);
        N1:=Filtered(N1,x->not x in nc);
    od;
    N2:=NormalSubgroups2(G2);
    nn2:=[];  Nnn2:=[];
    while N2<>[] do
        n:=N2[1];
        nn:=Unique(GeneratorsOfGroup(Normalizer(Group(Ng2),n)));
        Add(nn2,n);
        Add(Nnn2,nn);
        nc:=OuterConjugateSubgroups(Group(Ng2),n);
        N2:=Filtered(N2,x->not x in nc);
    od;
    ans:=[];
    for n1 in nn1 do
        epi1:=NaturalHomomorphismByNormalSubgroup(G1,n1);
        for n2 in nn2 do
            epi2:=NaturalHomomorphismByNormalSubgroup(G2,n2);
            iso:=IsomorphismGroups(Image(epi1),Image(epi2));
            if iso<>fail then
                aut:=AutomorphismGroup(Image(epi2));
                cc1:=List(Nnn1[Position(nn1,n1)],x->ConjugatorAutomorphism(G1,x));
                c1:=Group(List(cc1,x->InducedAutomorphism(iso,InducedAutomorphism(epi1,x))));
                cc2:=List(Nnn2[Position(nn2,n2)],x->ConjugatorAutomorphism(G2,x));
                c2:=Group(List(cc2,x->InducedAutomorphism(epi2,x)));
                d:=DoubleCosets2(aut,c1,c2);
                dr:=List(d,Representative);
                for f in dr do
                    Add(ans,SubdirectProduct(G1,G2,CompositionMapping(f,iso,epi1),epi2));
                od;
            fi;
        od;
    od;
    return ans;
end;

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;

p2m:= function(G)
    local gg,G1,G2;
    gg:=GeneratorsOfGroup(G);
    G1:=Group(List(gg,x->Image(Projection(G,1),x)));
    G2:=Group(List(gg,x->Image(Projection(G,2),x)));
    return DirectSumMatrixGroup([G1,G2]);
end;