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

InvariantFormsMatrix:= function(G)
    local d,gg,ut,bs,b,i,m,l,g,p;
    d:=Length(Identity(G));
    gg:=GeneratorsOfGroup(G);
    ut:=UnorderedTuples([1..d],2);
    bs:=[];
    for i in ut do
        b:=NullMat(d,d);
        b[i[1]][i[2]]:=1;
        b[i[2]][i[1]]:=1;
        Add(bs,b);
    od;
    if gg=[] then
        return bs;
    fi;
    m:=[];
    for b in bs do
        l:=[];
        for g in gg do
            p:=g*b*TransposedMat(g)-b;
            l:=Concatenation(l,List(ut,x->p[x[1]][x[2]]));
        od;
        Add(m,l);
    od;
    return List(NullspaceIntMat(m),x->x*bs);
end;

InvariantQuadraticForms:= function(G)
    local d,R,xx,gg,ut,bp,bm,b,i,bx,m,l,g,p;
    d:=Length(Identity(G));
    R:=PolynomialRing(Integers,d);
    xx:=IndeterminatesOfPolynomialRing(R);
    gg:=GeneratorsOfGroup(G);
    ut:=UnorderedTuples([1..d],2);
    bp:=List(ut,x->xx[x[1]]*xx[x[2]]);
    bm:=[];
    for i in ut do
        b:=NullMat(d,d);
        if i[1]=i[2] then
            b[i[1]][i[2]]:=1;
        else
            b[i[1]][i[2]]:=1/2;
            b[i[2]][i[1]]:=1/2;
        fi;
        Add(bm,b);
    od;
    if gg=[] then
        bx:=bp;
    else
        m:=[];
        for b in bm do
            l:=[];
            for g in gg do
                p:=g*b*TransposedMat(g)-b;
                l:=Concatenation(l,List(ut,x->p[x[1]][x[2]]*Length(Set(x))));
            od;
            Add(m,l);
        od;
        bx:=List(NullspaceIntMat(m),x->x*bp);
    fi;
    return bx;
end;

AutomorphismGroupOfMatrix:= function(M)
    local d,LM,M0,P,nb,m,vv,cb,ii,sp,gg,g;
    d:=Length(M);
    LM:=LLLReducedGramMat(M);
    M0:=LM.remainder;
    P:=LM.transformation;
    nb:=List([1..d],x->M0[x][x]);
    m:=Maximum(nb);
    vv:=ShortestVectors(M0,m);
    cb:=List(nb,x->Filtered(vv.vectors,y->y*M0*y=x));
    cb:=List(cb,x->Concatenation(List(x,y->[y,-y])));
    ii:=List([1..d],Zero);
    ii[1]:=1;
    sp:=1;
    gg:=[];
    while sp>0 do
        if ii[sp]>Length(cb[sp]) then
            sp:=sp-1;
            if sp>0 then
                ii[sp]:=ii[sp]+1;
            fi;
        elif ForAll([1..sp-1],x->cb[x][ii[x]]*M0*cb[sp][ii[sp]]=M0[x][sp]) then
            if sp<d then
                sp:=sp+1;
                ii[sp]:=1;
            else
                g:=List([1..d],x->cb[x][ii[x]])^P;
                if not g in Group(gg,IdentityMat(d)) then
                    Add(gg,g);
                fi;
                ii[sp]:=ii[sp]+1;
            fi;
        else
            ii[sp]:=ii[sp]+1;
        fi;
    od;
    return Group(gg);
end;

QuadraticFormToMatrix:= function(arg)
    local f,R,xx,d,M;
    f:=arg[1];
    if Length(arg)=1 then
        R:=DefaultRing(f);
        xx:=IndeterminatesOfPolynomialRing(R);
    elif IsRing(arg[2]) then
        R:=arg[2];
        xx:=IndeterminatesOfPolynomialRing(R);
    else
        xx:=arg[2];
    fi;
    d:=Length(xx);
    M:=List([1..d],i->List([1..d],j->Derivative(Derivative(f,xx[i]),xx[j])));
    M:=List([1..d],i->List([1..d],j->Value(M[i][j],xx,List(xx,x->0))));
    return M/2;
end;

MatrixToQuadraticForm:= function(arg)
    local M,d,R,xx;
    M:=arg[1];
    d:=Length(M);
    if Length(arg)=1 then
        R:=PolynomialRing(Integers,d);
        xx:=IndeterminatesOfPolynomialRing(R);
    elif IsRing(arg[2]) then
        R:=arg[2];
        xx:=IndeterminatesOfPolynomialRing(R);
    else
        xx:=arg[2];
    fi;
    return xx*M*xx;
end;

IsPositiveDefinite:= function(M)
    local d;
    d:=Length(M);
    return ForAll([1..d],x->Determinant(M{[1..x]}{[1..x]})>0);
end;