# The following code is used in
# "Rationality Problem for Algebraic Tori"
# by Akinari Hoshi and Aiichi Yamasaki
# Written by Aiichi Yamasaki
# The files in crystdat.zip are required.

RequirePackage("carat");

Read("cryst1.txt");
Read("cryst2.txt");
Read("cryst3.txt");
Read("cryst4.txt");
Read("cryst5.txt");
Read("cryst6.txt");

cryst:=[cryst1,cryst2,cryst3,cryst4,cryst5,cryst6];

Read("caratchpol.txt");

chpol:=[chpol1,chpol2,chpol3,chpol4,chpol5,chpol6];

Read("H1cryst.txt");

h1glnz:=[h1gl1z,h1gl2z,h1gl3z,h1gl4z,h1gl5z,h1gl6z];

Read("carat2crystcat.txt");

CharacteristicCyclotomicPolynomial:= function(m)
    local d,l,i;
    d:=Length(m);
    l:=[];
    for i in DivisorsInt(Order(m)) do
        l:=Concatenation(l,List([RankMat(m-E(i)*IdentityMat(d))+1..d],x->i));
    od;
    return l;
end;

CharPolyList:= function(g)
    local c;
    c:=List(ConjugacyClasses(g),Representative);
    return SortedList(List(c,CharacteristicCyclotomicPolynomial));
end;

CharPolySubgroupsList:= function(g)
    local h;
    h:=List(ConjugacyClassesSubgroups(g),Representative);
    return SortedList(List(h,CharPolyList));
end;

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

H1:= function(g)
    local m,gg,i,s,r;
    m:=[];
    gg:=GeneratorsOfGroup(g);
    if gg=[] then
        return [];
    else
        for i in gg do
            m:=Concatenation(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;

CaratQClass:= function(G)
    local d,ch,q;
    d:=Length(Identity(G));
    ch:=CharPolyList(G);
    q:=Position(chpol[d],ch);
    if d=6 and (q in chpol6dup) then
        ch:=CharPolySubgroupsList(G);
        q:=chpol6dup[Position(chpol6sub,ch)];
    fi;
    return [d,q];
end;

CaratZClass:= function(G)
    local cat,z,zz,h,hh,zg;
    cat:=CaratQClass(G);
    hh:=h1glnz[cat[1]][cat[2]];
    if Length(hh)=1 then
        return Concatenation(cat,[1]);
    else
        h:=[Hminus1(G),H0(G),H1(G)];
        zz:=Filtered([1..Length(hh)],x->hh[x]=h);
        if Length(zz)=1 then
            return Concatenation(cat,zz);
        else
            for z in zz do
                zg:=Group(cryst[cat[1]][cat[2]][z],IdentityMat(cat[1]));
                if RepresentativeAction(GL(cat[1],Integers),G,zg) in GL(cat[1],Integers) then
                    return Concatenation(cat,[z]);
                fi;
            od;
        fi;
    fi;
end;

CaratQClassNumber:= function(G)
    return CaratQClass(G);
end;

CaratZClassNumber:= function(G)
    return CaratZClass(G);
end;

CaratMatGroupZClass:= function(d,m,n)
    return Group(cryst[d][m][n],IdentityMat(d));
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;

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;

Carat2CrystCat:= function(l)
    local d,c;
    d:=l[1];
    c:=carat2crystcatQ[d][l[2]];
    if Length(l)=3 then
        c:=Concatenation(c,[carat2crystcatZ[d][l[2]][l[3]]]);
    fi;
    return c;
end;

CrystCat2Carat:= function(l)
    local d,n,c;
    d:=l[1];
    if Length(l)=3 then
        n:=Position(carat2crystcatQ[d],l);
        return [d,n];
    else
        n:=Position(carat2crystcatQ[d],[l[1],l[2],l[3]]);
        c:=[d,n];
        c[3]:=Position(carat2crystcatZ[d][n],l[4]);
        return c;
    fi;
end;

CrystCatQClass:= function(G)
    return Carat2CrystCat(CaratQClass(G));
end;

CrystCatZClass:= function(G)
    return Carat2CrystCat(CaratZClass(G));
end;

CrystCatQClassCatalog:= function(G)
    return Carat2CrystCat(CaratQClass(G));
end;

CrystCatZClassCatalog:= function(G)
    return Carat2CrystCat(CaratZClass(G));
end;

CrystCatQClassNumber:= function(G)
    return Carat2CrystCat(CaratQClass(G));
end;

CrystCatZClassNumber:= function(G)
    return Carat2CrystCat(CaratZClass(G));
end;

NrQClasses:= function(d)
    return Length(cryst[d]);
end;

NrZClasses:= function(d,i)
    return Length(cryst[d][i]);
end;

CaratNrQClasses:= function(d)
    return Length(cryst[d]);
end;
    
CaratNrZClasses:= function(d,i)
    return Length(cryst[d][i]);
end;