InstallGlobalFunction(FactorizationNParts, # decompose a positive integer d into product of n integers
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);
###################################################


InstallGlobalFunction(CrystGFullBasis,    # search for a G-fullbasis of a G-lattice 
                                          # for the given crystallographic G generated by 
        # the set of generators S
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;


##########
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" ]

];
###############

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 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) then 
##
 ctr:=Sum(B_delta)/2-vect; #test if one vertex of fundamental domain is origin
 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;


###
ctr:=0*vect; #test if center of fundamental domain is origin
 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);