# Compute the product of two sets in a group

ProdSets := function ( A, B )
    local AB, a, b;

    AB := [];
    for a in A do
        for b in B do
            AddSet(AB, a*b);
        od;
    od;

    return AB;
end;

# Compute a factorization of a group into susbets

FactorizationIntoSubsets:=function(G, H, K, a)
	local BuildKnuthMat, Knuth, b, aq, bq, iH, iK, RHT, LK, LKT, D, coset, Els, Dict, M, Combs, FirstKnuthMat, Solution, Way, lvl, PartSols, KnuthMat, KnuthArray, A, B, i, j, k, s, colnum, x, comb, h;
	BuildKnuthMat:=function(M, aq, bq, comb) #this function builds matrix of 0s and 1s for Knuth's algorithm X
		local Sets, CurSet, FirstKnuthMat, l, x;
		Sets:=[]; #assigning collection of sets
		for i in [1..Size(M[1])] do #filling the collection with sets
			CurSet:=[];
			for j in [1..aq] do
				Add(CurSet, M[comb[j]][i]);
			od;
			Add(Sets, CurSet);
		od;
		FirstKnuthMat:=NullMat(aq*bq+1, Size(Sets)); #assigning matrix of 0s and 1s with additional column for rows numeration, where the first index is column and the second one is row
		FirstKnuthMat[1]:=[1..Size(Sets)]; #numeration of rows
		for i in [1..Size(Sets)] do #filling the matrix with 1s
			for j in Sets[i] do
				FirstKnuthMat[j+1][i]:=1;
			od;
		od;
		for i in [1..Size(Sets)] do #removing rows obtained from lists with repeating elements that are not sets
			if Size(SSortedList(Sets[i]))<Size(Sets[i]) then
				x:=0;
				for j in FirstKnuthMat[1] do
					x:=x+1;
					if j=i then
						for colnum in [1..Size(FirstKnuthMat)] do
							Remove(FirstKnuthMat[colnum], x);
						od;
						break;
					fi;
				od;
			fi;
		od;
		return FirstKnuthMat;
	end;
	
	Knuth:=function(KnuthMat, Way, lvl, PartSols) #function that reduce matrix of 0s and 1s
		local ColSums, DetCol, c, RowNum, ColNums, RowNums;
		if Size(KnuthMat)=1 then #checking if matrix is empty and hence solution is obtained
			Solution:=1;
		else
			Solution:=0;
			ColSums:=[]; #assigning set of sums of columns elements
			for i in [2..Size(KnuthMat)] do
				AddSet(ColSums, Sum(KnuthMat[i]));
			od;
			if 0 in ColSums then #checking if there are columns with 0s only
				if lvl=1 then #if there are zero columns at level 1, there are no exact covers for given universe and collecion
					Solution:=-1;
				else
					Way[lvl-1]:=Way[lvl-1]+1; #moving to the next branch at level lvl-1
					for i in [lvl..Size(Way)] do
						Way[i]:=1;
					od;
					lvl:=1;
				fi;
			else
				for colnum in [2..Size(KnuthMat)] do #choosing column with the least number of 1s
					if Sum(KnuthMat[colnum])=ColSums[1] then
						DetCol:=KnuthMat[colnum];
						break;
					fi;
				od;
				c:=0;
				for i in [1..Size(DetCol)] do #choosing row according to the given way
					if DetCol[i]=1 then
						c:=c+1;
						if c=Way[lvl] then
							RowNum:=i;
							AddSet(PartSols,KnuthMat[1][i]); #including coumn of submatrix of matrix of double cosets as partial solution
							break;
						fi;
					fi;
				od;
				if c<Way[lvl] then #case when all appropriate rows are considered
					if lvl=1 then #if there are no more appropriate rows at level 1, no solutions for this matrix
						Solution:=-1;
					else
						Way[lvl-1]:=Way[lvl-1]+1; #moving to the next branch at level lvl-1
						for i in [lvl..Size(Way)] do
							Way[i]:=1;
						od;
						lvl:=1;
					fi;
				else
					ColNums:=[]; #picking numbers of columns that are to remove
					for i in [2..Size(KnuthMat)] do
						if KnuthMat[i][RowNum]=1 then
							AddSet(ColNums, i);
						fi;
					od;
					RowNums:=[]; #picking numbers of rows that are to remove
					for i in ColNums do
						for j in [1..Size(DetCol)] do
							if KnuthMat[i][j]=1 then
								AddSet(RowNums, j);
							fi;
						od;
					od;
					for i in RowNums do #removing rows
						for x in [1..Size(KnuthMat)] do
							Remove(KnuthMat[x], i);
						od;
						for j in [1..Size(RowNums)] do
							RowNums[j]:=RowNums[j]-1;
						od;
					od;
					for i in ColNums do #removing columns
						Remove(KnuthMat, i);
						for j in [1..Size(ColNums)] do
							ColNums[j]:=ColNums[j]-1;
						od;
					od;
					lvl:=lvl+1; #moving to the next level
					if Size(Way)<lvl then
						Add(Way, 1);
					fi;
				fi;
			fi;
		fi;
		return([KnuthMat, Way, lvl, Solution, PartSols]);
	end;

	b:=Order(G)/a;
	if not IsInt(b) then
		Print("Error, a does not divide order of G.");
		return;
	fi;
	aq:=a/Order(H);
	if not IsInt(aq) then
		Print("Error, order of H does not divide a.");
		return;
	fi;
	bq:=b/Order(K);
	if not IsInt(aq) then
		Print("Error, order of K does not divide b.");
		return;
	fi;
	iH:=Order(G)/Order(H);
	iK:=Order(G)/Order(K);
	RHT:=RightTransversal(G,H);
	LK:=DoubleCosets(G,TrivialSubgroup(G),K);
	LKT:=List([1..iK],x->Representative(LK[x]));
	D:=DoubleCosets(G,H,K);
	M:=NullMat(iH,iK);
	Els := Elements(G);;
    	Dict := NewDictionary(Els[1], true);;
    	k := 1;;
    	for coset in D do
       		for x in Elements(coset) do
         		AddDictionary(Dict, x, k);;
        	od;
        	k := k+1;;
    	od;
	for i in [1..iH] do
		for j in [1..iK] do
			M[i][j]:=LookupDictionary(Dict, RHT[i]*LKT[j]);
		od;
	od;
	Print("Matrix of double cosets was built.");
	Combs:=IteratorOfCombinations([1..iH], aq);
	for comb in Combs do
		FirstKnuthMat:=BuildKnuthMat(M, aq, bq, comb);
		Solution:=0;
		Way:=[1];
		lvl:=1;
		PartSols:=[];
		while Solution=0 do
			if lvl=1 then	
				KnuthMat:=StructuralCopy(FirstKnuthMat);
				PartSols:=[];
			fi;
			KnuthArray:=Knuth(KnuthMat, Way, lvl, PartSols);
			KnuthMat:=KnuthArray[1];
			Way:=KnuthArray[2];
			lvl:=KnuthArray[3];
			Solution:=KnuthArray[4];
			PartSols:=KnuthArray[5];
		od;
		if Solution=1 then
			break;
		fi;
	od;
	A:=[];
	B:=[];
	for i in comb do
		for h in H do
			AddSet(A,h*RHT[i]);
		od;
	od;
	for j in PartSols do
		for k in K do
			AddSet(B,LKT[j]*k);
		od;
	od;
	return [A, B];
end;