##########################################################################
#0
#F  FactorizationNParts
##
##  Decompose a positive integer d into product of n integers
##
##  Input: An pair of positive integers (d,n) 
##         
##  Output: A list of all possible ways to decompose d into 
##          product of n integers   
##
InstallGlobalFunction(FactorizationNParts, 
function(d,n)
local 
    t,faclst,x,w,y,j;
      
    if n=1 then 
        return [[d]];
    fi;
    
    t:=DivisorsInt(d);
    faclst:=[];
    for x in t do
        y:=FactorizationNParts(d/x,n-1);
        for j in y do
            Add(faclst,AddFirst(j,x));
        od;
    od;
    return faclst;
end);
################### end of FactorizationNParts ###########################

##########################################################################
#0
#F  CrystGFullBasis
##
##  Search for a G-fullbasis of a G-lattice
##  for the given crystallographic G generated by
##  the set of generators.
##
##  Input: a crystallographic group G 
##         
##  Output: If $G$ admits a parallelepiped fundamental region 
##          it outputs a $G$-full basis for a $G$-lattice $L$. 
##          If $G$ goes through the check for not admitting a 
##          parallelepiped fundamental region it outputs "false". 
##          Otherwise it outputs "fail".
##             
##
InstallGlobalFunction(CrystGFullBasis,                                              
function(arg)
local 
    G,P,Bt,d,n,vect,
    i,j,a,faclst,S,gens,L,c,
    B_delta,ctr,v,
    x,SbGrp,T,coef,SubgroupsOfAutCube;


# Create a list of subgroup of the automorphism group of an n-cube
    SubgroupsOfAutCube:=[

    [ "1", "C2", "C2 x C2", "C4", "D8" ],

    [ "1", "C2", "C3", "C2 x C2", "C4", "C6", "S3", "C2 x C2 x C2", 
    "D8", "C4 x C2", "A4", "D12", "C2 x D8", "C2 x A4", "S4", "C2 x S4" ],

    [ "1", "C2", "C3", "C2 x C2", "C4", "S3", "C6", "C2 x C2 x C2", "D8", 
    "C4 x C2", "C8", "Q8", "D12", "A4", "C6 x C2", "(C4 x C2) : C2", 
    "C2 x D8", "C2 x C2 x C2 x C2", "C4 x C4", "C8 : C2", "C4 : C4", 
    "D16", "QD16", "C4 x C2 x C2", "C2 x A4", "S4", "SL(2,3)", 
    "C2 x C2 x S3", "C2 x C2 x D8", "(C4 x C2 x C2) : C2", "C4 x D8",
    "((C4 x C2) : C2) : C2", "(C2 x D8) : C2", "(C2 x C2 x C2 x C2) : C2", 
    "(C8 : C2) : C2", "(C4 x C4) : C2", "C2 x S4", "C2 x C2 x A4", 
    "GL(2,3)", "((C4 x C4) : C2) : C2", "(((C4 x C2) : C2) : C2) : C2", 
    "D8 x D8", "((C8 : C2) : C2) : C2", "C2 x C2 x S4", 
    "((C2 x D8) : C2) : C3", 
    "(D8 x D8) : C2", "(((C2 x D8) : C2) : C3) : C2", 
    "((((C2 x D8) : C2) : C3) : C2) : C2" ]

    ];
    ####################end of Subgroup of Aut(cube)##########################

    gens:=GeneratorsOfGroup(arg[1]);
    G:=AffineCrystGroup(gens);
    SbGrp:=TranslationSubGroup(G);
    Bt:=G!.TranslationBasis;
    n:=G!.DimensionOfMatrixGroup-1;
    P:=PointGroup(G);

    if not StructureDescription(P) in SubgroupsOfAutCube[n-1] then
        return false;
    fi;

    if not GramianOfAverageScalarProductFromFiniteMatrixGroup(P)=
                                               IdentityMat(n) then
        return "Gramian matrix is not identity matrix";
    fi;

    if Length(arg)=1 then 
        L:=CrystCubicalTiling(n);
        Add(L,Bt,1);
        for i in [1..Length(L)] do
            B_delta:=CrystGFullBasis(G,[L[i],Sum(L[i])/2]);
            if IsList(B_delta) then 
                Add(B_delta,L[i]);   #Added October 2024
                return B_delta;
            fi;
        od;
        return fail;
    else #begin of length 2 input
        L:=arg[2][1];
        c:=arg[2][2];
        vect:=Sum(L)/2-c;
        S:=RightTransversal(G,SbGrp);
        d:=DivisorsInt(Order(P));
        i:=Length(d);

        while i>0 do
            faclst:=FactorizationNParts(d[i],n);
            for x in faclst do
                B_delta:=List([1..n],i->x[i]^-1*L[i]);
                if IsCrystSufficientLattice(B_delta,S,SbGrp) then 
                
 #test if one vertex of fundamental domain is origin
                    ctr:=Sum(B_delta)/2-vect; 
                    coef:=CombinationDisjointSets(x);
                    j:=1;
                    while j <= d[i] do
                        v:=ctr+coef[j]*B_delta;
                        if IsCrystSameOrbit(G,Bt,S,ctr,v)=false then 
                            break;
                        fi;
                        j:=j+1;
                    od;
                    if j=d[i]+1 then 
                        return [B_delta,ctr]; 
                    fi;

#test if center of fundamental domain is origin
                    ctr:=0*vect; 
                    j:=1;
                    while j <= d[i] do
                        v:=ctr+coef[j]*B_delta;
                        if IsCrystSameOrbit(G,Bt,S,ctr,v)=false then 
                            break;
                        fi;
                        j:=j+1;
                    od;
                    if j=d[i]+1 then 
                        return [B_delta,ctr];  
                    fi;
                fi;
            od;
            i:=i-1;
        
        od;
        return fail;
    fi;  #end of length 2 input

end);
################### end of CrystGFullBasis ###############################