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

Zminus1:= function(g)
    local m,gg,i,s;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        for i in gg do
            Append(m,i-Identity(g));
        od;
        s:=SmithNormalFormIntegerMatTransforms(m);
        m:=Inverse(s.coltrans);
        return List([1..s.rank],x->m[x]);
    fi;
end;

Bminus1:= function(g)
    local m,gg,i;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    for i in gg do
        Append(m,i-Identity(g));
    od;
    return LatticeBasis(m);
end;

Hminus1:= function(g)
    local m,gg,i,s,r;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        for i in gg do
            Append(m,i-Identity(g));
        od;
        s:=SmithNormalFormIntegerMat(m);
        r:=Rank(s);
        return List([1..r],x->s[x][x]);
    fi;
end;

Z0:= function(g)
    local m,s;
    m:=Sum(g);
    s:=SmithNormalFormIntegerMatTransforms(m);
    m:=Inverse(s.coltrans);
    return List([1..s.rank],x->m[x]);
end;

B0:= function(g)
    return LatticeBasis(Sum(g));
end;    

H0:= function(g)
    local m,s,r;
    m:=Sum(g);
    s:=SmithNormalFormIntegerMat(m);
    r:=Rank(s);
    return List([1..r],x->s[x][x]);
end;

Z1:= function(g)
    local m,gg,i,s;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        for i in gg do
            Append(m,TransposedMat(i)-Identity(g));
        od;
        m:=TransposedMat(m);
        s:=SmithNormalFormIntegerMatTransforms(m);
        m:=Inverse(s.coltrans);
        return List([1..s.rank],x->m[x]);
    fi;
end;

B1:= function(g)
    local m,gg,i;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    for i in gg do
        Append(m,TransposedMat(i)-Identity(g));
    od;
    m:=TransposedMat(m);
    return LatticeBasis(m);
end;

H1:= function(g)
    local m,gg,i,s,r;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        for i in gg do
            Append(m,TransposedMat(i)-Identity(g));
        od;
        m:=TransposedMat(m);
        s:=SmithNormalFormIntegerMat(m);
        r:=Rank(s);
        return List([1..r],x->s[x][x]);
    fi;
end;

Z2:= function(g)
    local gg,o,eg,j1,j2,j,l0,l1,l2,m0,m1,m2,m,s;
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return Identity(g);
    else
        o:=Order(g);
        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,o,o*Length(gg));
        m1:=BlockMatrix(l1,o,o*Length(gg));
        m2:=BlockMatrix(l2,o,o*Length(gg));
        m:=MatrixByBlockMatrix(m0-m1+m2);
        s:=SmithNormalFormIntegerMatTransforms(m);
        m:=Inverse(s.coltrans);
        return List([1..s.rank],x->m[x]);
    fi;
end;

B2:= function(g)
    local gg,o,eg,j1,j2,j,l0,l1,l2,m0,m1,m2,m;
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return Identity(g);
    else
        o:=Order(g);
        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,o,o*Length(gg));
        m1:=BlockMatrix(l1,o,o*Length(gg));
        m2:=BlockMatrix(l2,o,o*Length(gg));
        m:=MatrixByBlockMatrix(m0-m1+m2);
        return LatticeBasis(m);
    fi;
end;

H2:= function(g)
    local gg,o,eg,j1,j2,j,l0,l1,l2,m0,m1,m2,m,s,r;
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return List([1..Length(Identity(g))],x->1);
    else
        o:=Order(g);
        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,o,o*Length(gg));
        m1:=BlockMatrix(l1,o,o*Length(gg));
        m2:=BlockMatrix(l2,o,o*Length(gg));
        m:=MatrixByBlockMatrix(m0-m1+m2);
        s:=SmithNormalFormIntegerMat(m);
        r:=Rank(s);
        return List([1..r],x->s[x][x]);
    fi;
end;

Z3:= function(g)
    local gg,o,eg,j1,j2,j3,j,i1,i2,l0,l1,l2,l3,m0,m1,m2,m3,m,s;
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        o:=Order(g);
        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)*o+i2,j,Identity(g)]);
                    i1:=Position(eg,j1*j2);
                    Add(l1,[(i1-1)*o+i2,j,Identity(g)]);
                    i1:=Position(eg,j1); i2:=Position(eg,j2*j3);
                    Add(l2,[(i1-1)*o+i2,j,Identity(g)]);
                    i2:=Position(eg,j2);
                    Add(l3,[(i1-1)*o+i2,j,j3]);
                od;
            od;
        od;
        m0:=BlockMatrix(l0,o^2,o^2*Length(gg));
        m1:=BlockMatrix(l1,o^2,o^2*Length(gg));
        m2:=BlockMatrix(l2,o^2,o^2*Length(gg));
        m3:=BlockMatrix(l3,o^2,o^2*Length(gg));
        m:=MatrixByBlockMatrix(m0-m1+m2-m3);
        s:=SmithNormalFormIntegerMatTransforms(m);
        m:=Inverse(s.coltrans);
        return List([1..s.rank],x->m[x]);
    fi;
end;

B3:= function(g)
    local gg,o,eg,j1,j2,j3,j,i1,i2,l0,l1,l2,l3,m0,m1,m2,m3,m;
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        o:=Order(g);
        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)*o+i2,j,Identity(g)]);
                    i1:=Position(eg,j1*j2);
                    Add(l1,[(i1-1)*o+i2,j,Identity(g)]);
                    i1:=Position(eg,j1); i2:=Position(eg,j2*j3);
                    Add(l2,[(i1-1)*o+i2,j,Identity(g)]);
                    i2:=Position(eg,j2);
                    Add(l3,[(i1-1)*o+i2,j,j3]);
                od;
            od;
        od;
        m0:=BlockMatrix(l0,o^2,o^2*Length(gg));
        m1:=BlockMatrix(l1,o^2,o^2*Length(gg));
        m2:=BlockMatrix(l2,o^2,o^2*Length(gg));
        m3:=BlockMatrix(l3,o^2,o^2*Length(gg));
        m:=MatrixByBlockMatrix(m0-m1+m2-m3);
        return LatticeBasis(m);
    fi;
end;

H3:= function(g)
    local gg,o,eg,j1,j2,j3,j,i1,i2,l0,l1,l2,l3,m0,m1,m2,m3,m,s,r;
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        o:=Order(g);
        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)*o+i2,j,Identity(g)]);
                    i1:=Position(eg,j1*j2);
                    Add(l1,[(i1-1)*o+i2,j,Identity(g)]);
                    i1:=Position(eg,j1); i2:=Position(eg,j2*j3);
                    Add(l2,[(i1-1)*o+i2,j,Identity(g)]);
                    i2:=Position(eg,j2);
                    Add(l3,[(i1-1)*o+i2,j,j3]);
                od;
            od;
        od;
        m0:=BlockMatrix(l0,o^2,o^2*Length(gg));
        m1:=BlockMatrix(l1,o^2,o^2*Length(gg));
        m2:=BlockMatrix(l2,o^2,o^2*Length(gg));
        m3:=BlockMatrix(l3,o^2,o^2*Length(gg));
        m:=MatrixByBlockMatrix(m0-m1+m2-m3);
        s:=SmithNormalFormIntegerMat(m);
        r:=Rank(s);
        return List([1..r],x->s[x][x]);
    fi;
end;

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

IsFlabby:= function(g)
    local h;
    h:=List(ConjugacyClassesSubgroups2(g),Representative);
    return ForAll(h,x->Product(Hminus1(x))=1);
end;

IsCoflabby:= function(g)
    local h;
    h:=List(ConjugacyClassesSubgroups2(g),Representative);
    return ForAll(h,x->Product(H1(x))=1);
end;