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

BlockList:= function(L,n)
    local l,m;
    l:=Length(L);
    m:=l/n;
    return List([1..n],x->L{[(x-1)*m+1..x*m]});
end;

IsBicyclic:= function(g)
    local f3,p;
    if not IsAbelian(g) then
        return false;
    fi;
    f3:=Filtered(Collected(Factors(Order(g))),x->x[2]>2);
    for p in f3 do
        if Length(AbelianInvariants(SylowSubgroup(g,p[1])))>2 then
            return false;
        fi;
    od;
    return true;
end;

FindGenFiniteAbelian:= function(g)
    local e,a,ga,iso;
    e:=AbelianInvariants(g);
    if Length(e)>1 then
        e:=SmithNormalFormIntegerMat(DiagonalMat(e));
        e:=List([1..Length(e)],x->e[x][x]);
        e:=Filtered(e,x->x>1);
    fi;
    a:=AbelianGroup(e);
    ga:=GeneratorsOfGroup(a);
    iso:=IsomorphismGroups(a,g);
    return List(ga,x->Image(iso,x));
end;

EltFiniteAbelian:= function(g,c)
    local gg,F,gF,hom,cF,e;
    gg:=GeneratorsOfGroup(g);
    F:=FreeGroup(Length(gg));
    gF:=GeneratorsOfGroup(F);
    hom:=GroupHomomorphismByImages(F,g,gF,gg);
    cF:=PreImagesRepresentative(hom,c);
    e:=List(gF,x->ExponentSumWord(cF,x));
    return e;
end;

CheckSNF:= function(m)
    local s,r;
    s:=SmithNormalFormIntegerMat(m);
    r:=Rank(s);
    if r=0 then
        return 0;
    else
        return s[r][r];
    fi;
end;

AdjustListLength:= function(l,n)
    if Length(l)<n then
        return Concatenation(l,List([Length(l)+1..n],Zero));
    else
        return l{[1..n]};
    fi;
end;

ChangeHn:= function(f,Bn,zero,s)
    local l;
    l:=f{zero}*s.coltrans;
    l:=AdjustListLength(l,Length(Bn));
    return f-l*s.rowtrans*Bn;
end;

Z2value:= function(z,eg,gg)
    local og,zmat,i1,i2,i3,i,j2,j3,s,s1,s2;
    og:=Length(eg);
    zmat:=NullMat(og,og);
    for j2 in [1..Length(gg)] do
        i2:=Position(eg,gg[j2]);
        zmat{[1..og]}{[i2]}:=z{[1..og]}{[j2]};
    od;
    s:=gg;
    s1:=gg;
    s2:=[];
    repeat
        for j2 in s1 do
            i2:=Position(eg,j2);
            for j3 in gg do
                if not j2*j3 in Concatenation(s,s2) then
                    Add(s2,j2*j3);
                    i3:=Position(eg,j3);
                    i:=Position(eg,j2*j3);
                    for i1 in [1..og] do
                        zmat[i1][i]:=zmat[i1][i2]*j3
                            +zmat[Position(eg,eg[i1]*j2)][i3]-zmat[i2][i3];
                    od;
                fi;
            od;
        od;
        s1:=s2;
        s:=Union(s,s1);
        s2:=[];
    until s1=[];
    return zmat;
end;

Z3value:= function(z,eg,gg)
    local og,zmat,i1,i2,i3,i4,i,j3,j4,s,s1,s2;
    og:=Length(eg);
    zmat:=List([1..og],x->NullMat(og,og));
    for j3 in [1..Length(gg)] do
        i3:=Position(eg,gg[j3]);
        zmat{[1..og]}{[1..og]}{[i3]}:=z{[1..og]}{[1..og]}{[j3]};
    od;
    s:=gg;
    s1:=gg;
    s2:=[];
    repeat
        for j3 in s1 do
            i3:=Position(eg,j3);
            for j4 in gg do
                if not j3*j4 in Concatenation(s,s2) then
                    Add(s2,j3*j4);
                    i4:=Position(eg,j4);
                    i:=Position(eg,j3*j4);
                    for i1 in [1..og] do
                        for i2 in [1..og] do
                            zmat[i1][i2][i]:=zmat[i1][i2][i3]*j4
                                +zmat[i1][Position(eg,eg[i2]*j3)][i4]
                                -zmat[Position(eg,eg[i1]*eg[i2])][i3][i4]
                                +zmat[i2][i3][i4];
                        od;
                    od;
                fi;
            od;
        od;
        s1:=s2;
        s:=Union(s,s1);
        s2:=[];
    until s1=[];
    return zmat;
end;

ResH2:= function(g,h)
    local d,gg,og,eg,gh,oh,eh,j1,j2,j,i1,i2,i,l0,l1,l2,m0,m1,m2,m,
        zero,sg,sh,h2g,h2gg,h2h,h2hg,z,zmat,s,res,ress,res0,res0s;
    d:=Length(Identity(g));
    gg:=GeneratorsOfGroup(g);
    og:=Order(g);
    if gg=[] or og=1 then
        return rec(H2G:=[],H2Ggen:=[],H2H:=[],H2Hgen:=[],ResMat:=[]);
    fi;
    eg:=Elements(g);
    l0:=[]; l1:=[]; l2:=[];
    j:=0;
    for j1 in eg do
        for j2 in gg do
            j:=j+1;
            Add(l0,[Position(eg,j2),j,Identity(g)]);
            Add(l1,[Position(eg,j1*j2),j,Identity(g)]);
            Add(l2,[Position(eg,j1),j,j2]);
        od;
    od;
    m0:=BlockMatrix(l0,og,og*Length(gg));
    m1:=BlockMatrix(l1,og,og*Length(gg));
    m2:=BlockMatrix(l2,og,og*Length(gg));
    m:=MatrixByBlockMatrix(m0-m1+m2);
    zero:=Flat(List([0..Order(gg[1])-2],
        x->List([1..d],y->(Position(eg,gg[1]^x)-1)*Length(gg)*d+y)));
    m:=NullspaceIntMat(m{[1..Length(m)]}{zero})*m;
    m:=LatticeBasis(m);
    z:=1;
    while m<>[] and z<=og*Length(gg) do
        zero:=[(z-1)*d+1..z*d];
        l0:=m{[1..Length(m)]}{zero};
        if CheckSNF(l0)=1 then
            m:=NullspaceIntMat(l0)*m;
        fi;
        z:=z+1;
    od;
    if m=[] then
        h2g:=[];
        h2gg:=[];
    else
        sg:=SmithNormalFormIntegerMatTransforms(m);
        h2g:=Filtered(List([1..sg.rank],x->sg.normal[x][x]),y->y>1);
        m:=Inverse(sg.coltrans);
        h2gg:=List([sg.rank-Length(h2g)+1..sg.rank],x->m[x]);
        h2gg:=List(h2gg,x->BlockList(x,og));
        h2gg:=List(h2gg,x->List(x,y->BlockList(y,Length(gg))));
    fi;
    gh:=GeneratorsOfGroup(h);
    oh:=Order(h);
    if gh=[] or oh=1 then
        return rec(H2G:=h2g,H2Ggen:=h2gg,H2H:=[],H2Hgen:=[],
            ResMat:=List(h2g,x->[]));
    fi;
    eh:=Elements(h);
    l0:=[]; l1:=[]; l2:=[];
    j:=0;
    for j1 in eh do
        for j2 in gh do
            j:=j+1;
            Add(l0,[Position(eh,j2),j,Identity(h)]);
            Add(l1,[Position(eh,j1*j2),j,Identity(h)]);
            Add(l2,[Position(eh,j1),j,j2]);
        od;
    od;
    m0:=BlockMatrix(l0,oh,oh*Length(gh));
    m1:=BlockMatrix(l1,oh,oh*Length(gh));
    m2:=BlockMatrix(l2,oh,oh*Length(gh));
    m:=MatrixByBlockMatrix(m0-m1+m2);
    ress:=[];
    for z in h2gg do
        zmat:=Z2value(z,eg,gg);
        res:=List(eh,x->List(gh,y->zmat[Position(eg,x)][Position(eg,y)]));
        Add(ress,Flat(res));
    od;
    res0s:=[];
    for res in ress do
        res0:=res;
        i:=Position(eh,gh[1]);
        for i1 in [1..Order(gh[1])-2] do
            l0:=List([1..d],x->res0[(i-1)*Length(gh)*d+x]);
            i:=Position(eh,gh[1]^(i1+1));
            res0:=res0+l0*m{[(i-1)*d+1..i*d]};
        od;
        Add(res0s,res0);
    od;
    zero:=Flat(List([0..Order(gh[1])-2],
        x->List([1..d],y->(Position(eh,gh[1]^x)-1)*Length(gh)*d+y)));
    m:=NullspaceIntMat(m{[1..Length(m)]}{zero})*m;
    m:=LatticeBasis(m);
    z:=1;
    while m<>[] and z<=oh*Length(gh) do
        zero:=[(z-1)*d+1..z*d];
        l0:=m{[1..Length(m)]}{zero};
        if CheckSNF(l0)=1 then
            s:=SmithNormalFormIntegerMatTransforms(l0);
            res0s:=List(res0s,x->ChangeHn(x,m,zero,s));
            m:=NullspaceIntMat(l0)*m;
        fi;
        z:=z+1;
    od;
    if m=[] then
        return rec(H2G:=h2g,H2Ggen:=h2gg,H2H:=[],H2Hgen:=[],
            ResMat:=List(h2g,x->[]));
    fi;
    sh:=SmithNormalFormIntegerMatTransforms(m);
    h2h:=Filtered(List([1..sh.rank],x->sh.normal[x][x]),y->y>1);
    m:=Inverse(sh.coltrans);
    h2hg:=List([sh.rank-Length(h2h)+1..sh.rank],x->m[x]);
    h2hg:=List(h2hg,x->BlockList(x,oh));
    h2hg:=List(h2hg,x->List(x,y->BlockList(y,Length(gh))));
    ress:=[];
    for res0 in res0s do
        res:=res0*sh.coltrans;
        res:=List([sh.rank-Length(h2h)+1..sh.rank],
            x->res[x] mod sh.normal[x][x]);
        Add(ress,res);
    od;
    return rec(H2G:=h2g,H2Ggen:=h2gg,H2H:=h2h,H2Hgen:=h2hg,ResMat:=ress);
end;

H2nrM:= function(g)
    local d,gg,og,eg,h,gh,oh,eh,j1,j2,j,i1,i2,l0,l1,l2,m0,m1,m2,m,
        zero,sg,sh,h2g,h2gg,h2h,z,zmat,zmats,s,res,ress,
        Hg,Hgg,K,ga,iso,Hh,Hhg,hom,Kg;
    d:=Length(Identity(g));
    gg:=GeneratorsOfGroup(g);
    og:=Order(g);
    if gg=[] or og=1 then
        return rec(H2G:=[],H2Ggen:=[],H2nrM:=[],H2nrMgen:=[]);
    fi;
    eg:=Elements(g);
    l0:=[]; l1:=[]; l2:=[];
    j:=0;
    for j1 in eg do
        for j2 in gg do
            j:=j+1;
            Add(l0,[Position(eg,j2),j,Identity(g)]);
            Add(l1,[Position(eg,j1*j2),j,Identity(g)]);
            Add(l2,[Position(eg,j1),j,j2]);
        od;
    od;
    m0:=BlockMatrix(l0,og,og*Length(gg));
    m1:=BlockMatrix(l1,og,og*Length(gg));
    m2:=BlockMatrix(l2,og,og*Length(gg));
    m:=MatrixByBlockMatrix(m0-m1+m2);
    zero:=Flat(List([0..Order(gg[1])-2],
        x->List([1..d],y->(Position(eg,gg[1]^x)-1)*Length(gg)*d+y)));
    m:=NullspaceIntMat(m{[1..Length(m)]}{zero})*m;
    m:=LatticeBasis(m);
    z:=1;
    while m<>[] and z<=og*Length(gg) do
        zero:=[(z-1)*d+1..z*d];
        l0:=m{[1..Length(m)]}{zero};
        if CheckSNF(l0)=1 then
            m:=NullspaceIntMat(l0)*m;
        fi;
        z:=z+1;
    od;
    if m=[] then
        return rec(H2G:=[],H2Ggen:=[],H2nrM:=[],H2nrMgen:=[]);
    fi;
    sg:=SmithNormalFormIntegerMatTransforms(m);
    h2g:=Filtered(List([1..sg.rank],x->sg.normal[x][x]),y->y>1);
    m:=Inverse(sg.coltrans);
    h2gg:=List([sg.rank-Length(h2g)+1..sg.rank],x->m[x]);
    h2gg:=List(h2gg,x->BlockList(x,og));
    h2gg:=List(h2gg,x->List(x,y->BlockList(y,Length(gg))));
    zmats:=List(h2gg,x->Z2value(x,eg,gg));
    Hg:=AbelianGroup(h2g);
    Hgg:=GeneratorsOfGroup(Hg);
    K:=Hg;
    if ValueOption("fromperm")=true or ValueOption("FromPerm")=true then
        iso:=IsomorphismPermGroup(g);
        ga:=List(ConjugacyClassesSubgroups(Range(iso)),Representative);
        ga:=Filtered(ga,IsBicyclic);
        ga:=Concatenation(List(ga,x->ConjugateSubgroups(Image(iso),x)));
        ga:=List(ga,x->PreImage(iso,x));
    else
        ga:=List(ConjugacyClassesSubgroups(g),Representative);
        ga:=Filtered(ga,IsBicyclic);
        ga:=Concatenation(List(ga,x->ConjugateSubgroups(g,x)));
    fi;
    ga:=Filtered(ga,x->Order(x)>1);
    for h in ga do
        gh:=FindGenFiniteAbelian(h);
        oh:=Order(h);
        eh:=Elements(h);
        l0:=[]; l1:=[]; l2:=[];
        j:=0;
        for j1 in eh do
            for j2 in gh do
                j:=j+1;
                Add(l0,[Position(eh,j2),j,Identity(h)]);
                Add(l1,[Position(eh,j1*j2),j,Identity(h)]);
                Add(l2,[Position(eh,j1),j,j2]);
            od;
        od;
        m0:=BlockMatrix(l0,oh,oh*Length(gh));
        m1:=BlockMatrix(l1,oh,oh*Length(gh));
        m2:=BlockMatrix(l2,oh,oh*Length(gh));
        m:=MatrixByBlockMatrix(m0-m1+m2);
        sh:=SmithNormalFormIntegerMatTransforms(m);
        h2h:=Filtered(List([1..sh.rank],x->sh.normal[x][x]),y->y>1);
        if h2h<>[] then
            Hh:=AbelianGroup(h2h);
            Hhg:=GeneratorsOfGroup(Hh);
            ress:=[];
            for zmat in zmats do
                res:=List(eh,x->List(gh,
                    y->zmat[Position(eg,x)][Position(eg,y)]));
                res:=Flat(res)*sh.coltrans;
                res:=Product([1..Length(h2h)],
                    x->Hhg[x]^res[sh.rank-Length(h2h)+x]);
                Add(ress,res);
            od;
            hom:=GroupHomomorphismByImages(Hg,Hh,Hgg,ress);
            K:=Intersection(K,Kernel(hom));
            if Order(K)=1 then
                return rec(H2G:=h2g,H2Ggen:=h2gg,H2nrM:=[],H2nrMgen:=[]);
            fi;
        fi;
    od;
    Kg:=FindGenFiniteAbelian(K);
    return rec(H2G:=h2g,H2Ggen:=h2gg,H2nrM:=List(Kg,Order),
        H2nrMgen:=List(Kg,x->EltFiniteAbelian(Hg,x)));
end;

ResH3:= function(g,h)
    local d,gg,og,eg,gh,oh,eh,j1,j2,j3,j,i1,i2,i3,
        l0,l1,l2,l3,m0,m1,m2,m3,m,
        zero,sg,sh,h3g,h3gg,h3h,h3hg,z,zmat,s,res,ress;
    d:=Length(Identity(g));
    gg:=GeneratorsOfGroup(g);
    og:=Order(g);
    if gg=[] or og=1 then
        return rec(H3G:=[],H3Ggen:=[],H3H:=[],H3Hgen:=[],ResMat:=[]);
    fi;
    eg:=Elements(g);
    l0:=[]; l1:=[]; l2:=[]; l3:=[];
    j:=0;
    for j1 in eg do
        for j2 in eg do
            for j3 in gg do
                j:=j+1;
                i1:=Position(eg,j2); i2:=Position(eg,j3);
                Add(l0,[(i1-1)*og+i2,j,Identity(g)]);
                i1:=Position(eg,j1*j2);
                Add(l1,[(i1-1)*og+i2,j,Identity(g)]);
                i1:=Position(eg,j1); i2:=Position(eg,j2*j3);
                Add(l2,[(i1-1)*og+i2,j,Identity(g)]);
                i2:=Position(eg,j2);
                Add(l3,[(i1-1)*og+i2,j,j3]);
            od;
        od;
    od;
    m0:=BlockMatrix(l0,og^2,og^2*Length(gg));
    m1:=BlockMatrix(l1,og^2,og^2*Length(gg));
    m2:=BlockMatrix(l2,og^2,og^2*Length(gg));
    m3:=BlockMatrix(l3,og^2,og^2*Length(gg));
    m:=MatrixByBlockMatrix(m0-m1+m2-m3);
    z:=Position(eg,Identity(g));
    zero:=Union(List([1..og],x->[(x-1)*og+z,(z-1)*og+x]));
    zero:=Flat(List(zero,x->[(x-1)*Length(gg)*d+1..x*Length(gg)*d]));
    m:=NullspaceIntMat(m{[1..Length(m)]}{zero})*m;
    m:=LatticeBasis(m);
    z:=1;
    while m<>[] and z<=og^2*Length(gg) do
        zero:=[(z-1)*d+1..z*d];
        l0:=m{[1..Length(m)]}{zero};
        if CheckSNF(l0)=1 then
            m:=NullspaceIntMat(l0)*m;
        fi;
        z:=z+1;
    od;
    if m=[] then
        h3g:=[];
        h3gg:=[];
    else
        sg:=SmithNormalFormIntegerMatTransforms(m);
        h3g:=Filtered(List([1..sg.rank],x->sg.normal[x][x]),y->y>1);
        m:=Inverse(sg.coltrans);
        h3gg:=List([sg.rank-Length(h3g)+1..sg.rank],x->m[x]);
        h3gg:=List(h3gg,x->BlockList(x,og));
        h3gg:=List(h3gg,x->List(x,y->BlockList(y,og)));
        h3gg:=List(h3gg,x->List(x,y->List(y,z->BlockList(z,Length(gg)))));
    fi;
    gh:=GeneratorsOfGroup(h);
    oh:=Order(h);
    if gh=[] or oh=1 then
        return rec(H3G:=h3g,H3Ggen:=h3gg,H3H:=[],H3Hgen:=[],
            ResMat:=List(h3g,x->[]));
    fi;
    eh:=Elements(h);
    l0:=[]; l1:=[]; l2:=[]; l3:=[];
    j:=0;
    for j1 in eh do
        for j2 in eh do
            for j3 in gh do
                j:=j+1;
                i1:=Position(eh,j2); i2:=Position(eh,j3);
                Add(l0,[(i1-1)*oh+i2,j,Identity(h)]);
                i1:=Position(eh,j1*j2);
                Add(l1,[(i1-1)*oh+i2,j,Identity(h)]);
                i1:=Position(eh,j1); i2:=Position(eh,j2*j3);
                Add(l2,[(i1-1)*oh+i2,j,Identity(h)]);
                i2:=Position(eh,j2);
                Add(l3,[(i1-1)*oh+i2,j,j3]);
            od;
        od;
    od;
    m0:=BlockMatrix(l0,oh^2,oh^2*Length(gh));
    m1:=BlockMatrix(l1,oh^2,oh^2*Length(gh));
    m2:=BlockMatrix(l2,oh^2,oh^2*Length(gh));
    m3:=BlockMatrix(l3,oh^2,oh^2*Length(gh));
    m:=MatrixByBlockMatrix(m0-m1+m2-m3);
    z:=Position(eh,Identity(h));
    zero:=Union(List([1..oh],x->[(x-1)*oh+z,(z-1)*oh+x]));
    zero:=Flat(List(zero,x->[(x-1)*Length(gg)*d+1..x*Length(gg)*d]));
    m:=NullspaceIntMat(m{[1..Length(m)]}{zero})*m;
    m:=LatticeBasis(m);
    ress:=[];
    for z in h3gg do
        zmat:=Z3value(z,eg,gg);
        res:=List(eh,x->List(eh,y->List(gh,
            z->zmat[Position(eg,x)][Position(eg,y)][Position(eg,z)])));
        Add(ress,Flat(res));
    od;
    z:=1;
    while m<>[] and z<=oh^2*Length(gh) do
        zero:=[(z-1)*d+1..z*d];
        l0:=m{[1..Length(m)]}{zero};
        if CheckSNF(l0)=1 then
            s:=SmithNormalFormIntegerMatTransforms(l0);
            ress:=List(ress,x->ChangeHn(x,m,zero,s));
            m:=NullspaceIntMat(l0)*m;
        fi;
        z:=z+1;
    od;
    if m=[] then
        return rec(H3G1:=h3g,H3G1gen:=h3gg,H3G2:=[],H3G2gen:=[],
            ResMat:=List(h3g,x->[]));
    fi;
    sh:=SmithNormalFormIntegerMatTransforms(m);
    h3h:=Filtered(List([1..sh.rank],x->sh.normal[x][x]),y->y>1);
    m:=Inverse(sh.coltrans);
    h3hg:=List([sh.rank-Length(h3h)+1..sh.rank],x->m[x]);
    h3hg:=List(h3hg,x->BlockList(x,oh));
    h3hg:=List(h3hg,x->List(x,y->BlockList(y,oh)));
    h3hg:=List(h3hg,x->List(x,y->List(y,z->BlockList(z,Length(gh)))));
    ress:=List(ress,x->x*sh.coltrans);
    ress:=List(ress,x->List([sh.rank-Length(h3h)+1..sh.rank],
        y->x[y] mod sh.normal[y][y]));
    return rec(H3G:=h3g,H3Ggen:=h3gg,H3H:=h3h,H3Hgen:=h3hg,ResMat:=ress);
end;