(*{{2,3,4,9},{8,16}}
For this file, see the comments added to music16vs24839, the same apply here.

*)

Clear[x];
Expand[Cyclotomic[2,x]Cyclotomic[4,x]]

(*This section alone is the "fill-out procedure. You need to specify a value of "n" and a set 
"B" at the beginning and then run it. It gives you ALL the tiling complements of B in the group Z_n. Try to play with it in different groups. Z_108, Z_120, etc. It runs remarkably fast, as long as the number of tiling complements is not more than a few thousand. It takes longer if it is something like 100000, and takes 2-3 days if it is 1000000. Please use it at your wish. For example it is funny to give n=1000 and a random set B of length 20, and see how quickly it decides whether there are any tilining complements of B. Of course, almost surely B will have no tiling complements, and this procedure will run within 1-2 minutes. 
 *)

Dynamic[{i,g, mark}]
Clear[n,B,R,Tile,A,U,mark,T,g];

n=144;
B={0,1,2,3};
B=B⋃(B+16)⋃(B+32)⋃(B+48)⋃(B+64)⋃(B+80)⋃(B+96)⋃(B+112)⋃(B+128)

Tile=Array[0&,1000000];
R=Array[0&,n];
A=Array[200&,{n/Length[B],n,Length[B]}];
U=Array[0&,n/Length[B]];
U[[1]]=B;
mark=Array[1&,n/Length[B]+1];
T=Array[{0}&,n/Length[B]];
g=1;

For[i=1, i>0,
If[Length[U[[i]]]==n,  Tile[[g]]=T[[i]]⋃{};g=g+1;T[[i]]=T[[i-1]];mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];i=i-1,
For[j=1, j≤n,k=100;
If[Not[MemberQ[U[[i]],Mod[j,n]]], k=0;
For[m=1, m≤Length[B],
If[U[[i]]⋂ Mod[(j-B[[m]])+B,n]=={}, k=k+1; A[[i,j,m]]=Mod[j-B[[m]],n]];
m++]];
R[[j]]=k;
j++];
min=Min[R]; 
If[(Not[min==0]),
t=Catch[Do[If[R[[ii]]==min, Throw[ii]],{ii,n}]];
If[mark[[i]]≤Length[B],
For[b=mark[[i]], b<=Length[B],mark[[i]]=b+1;z=0;
If[Not[A[[i,t,b]]==200], U[[i+1]]=U[[i]]⋃Mod[A[[i,t,b]]+B,n]; T[[i+1]]=T[[i]]⋃{A[[i,t,b]]}; b=Length[B]+1;z=1;i=i+1;];
b++];
If[z==0, mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];i=i-1];
, mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];T[[i]]=T[[i-1]];i=i-1];
, mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];T[[i]]=T[[i-1]];i=i-1];
];
];

Tile=Drop[Tile⋃{},1];
Length[Tile]
Tile;

(* In this step we are going to throw away the periodic sets from Tile. And also, we want to keep only one copy of equivalent tiles (a tile is equivalent to its own translated coppies). We only want to keep the ones in normal form. The normal form i chose here is that the longest interval is thrown in at the beginning.

*)

For[j=1, j≤Length[Tile], P=Tile[[j]];
For[i=2, i<=Length[P], 
If[MemberQ[Tile,Mod[P-P[[i]],n]⋃{}], Tile[[j]]=0;i=Length[P]+1];
i++];
j++]; 
Tile=Drop[Tile⋃{},1];
Length[Tile]
Tile

(* This step is a seemingly long code, but it's task is simple. You must give a value "n" and a list of sets "T2" at the beginning as input. What it does: it puts the sets in T2 into equvalence classes according to which cyclotomic polynomials divide them. In this case n=144  and T2=Tile from the previous step.  
*)

n=144;

T2=Tile;
 
Dynamic[{j,v,vvv,vvvv}]


Clear[H];

Length[T2]

H=Array[0&,{Length[T2],n}];
MF=Array[0&,100];
MT=Array[0&,100];

For[j=1, j≤Length[T2], 
For[k=1, k≤n,
If[MemberQ[T2[[j]],Mod[k,n]],H[[j,k]]=1];
k++];
j++];

MF[[1]]=H[[1]];
MT[[1]]=T2[[1]];

cont=1;
v=1;
For[j=2, j≤Length[T2], cont=1;
For[c=1,c≤v,
For[k=1, k≤n,
If[(Not[Chop[Fourier[MF[[c]]]][[k]]==Chop[Fourier[H[[j]]]][[k]]])∧(Chop[Fourier[MF[[c]]]][[k]]==0∨Chop[Fourier[H[[j]]]][[k]]==0), 
k=n+1;cont=cont+1];
k++];
If[cont≤c, c=v+1];
c++];
If[cont==v+1, MF[[v+1]]=H[[j]];MT[[v+1]]=T2[[j]]; v=v+1];
j++]; 


MT=Drop[MT⋃{},1];
MF=Drop[MF⋃{},1];
Length[MF]



MT

 

(* Next, we see what those Fourier-zero sets are. *)

J=MF;
F1=Array[1&,{Length[J],n}];

For[i=1,i<=Length[J],
For[j=1,j≤n,
If[Chop[Fourier[J[[i]]]][[j]]==0, F1[[i,j]]=0];
j++];
i++];

F1;


FF1=Array[0&,{Length[F1],n}];

For[i=1,i≤Length[F1],
For[j=1,j≤n,
If[F1[[i,j]]==0, FF1[[i,j]]=j];
j++];
i++];

For[j=1,j≤Length[FF1],
FF1[[j]]=Drop[FF1[[j]]⋃{},1];
j++];
FF1;


FF=FF1-1;

For[m=1,m≤Length[FF1],
For[j=1,j≤Length[FF1[[m]]],
FF[[m,j]]=GCD[FF[[m,j]],n];
j++];
m++];
FF;
For[j=1,j≤Length[FF],
FF[[j]]=FF[[j]]⋃{};
j++];
FF;

FT=n/FF
CFT=Array[0&,Length[FT]];
For[j=1,j≤Length[FT],
CFT[[j]]=Complement[Divisors[n],FT[[j]]];
j++];
CFT

Dynamic[{jj,j,i,g, mark, rr}]
Clear[n,J,B,R,Tile,A,U,mark,T,g,x];

n=144;

J={{0,40,76,116},{0,44,84,104},{0,52,76,120},{0,52,88,108},{0,60,84,120},{0,84,108,120}};

rr=1;
gcd=Array[0&,{Length[J],Length[J[[1]]]}];
For[kk=1,kk≤Length[J],
For[mm=2,mm≤Length[J[[kk]]],
gcd[[kk,mm]]=GCD[gcd[[kk,mm-1]],J[[kk,mm]]];
mm++];
kk++];
For[kk=1,kk≤Length[J], gcd[[kk,Length[J[[1]]]]]=GCD[gcd[[kk,Length[J[[1]]]]],n];
kk++];

maxperiod=Array[0&,{Length[J],10000}];

For[jj=1, jj≤Length[J],Clear[B,R,Tile,A,U,mark,T,g,x];n=144;n=n/gcd[[jj, Length[J[[1]]]]];B=J[[jj]]/gcd[[jj,Length[J[[1]]]]];

Tile=Array[0&,100000];
R=Array[0&,n];
A=Array[200&,{n/Length[B],n,Length[B]}];
U=Array[0&,n/Length[B]];
U[[1]]=B;
mark=Array[1&,n/Length[B]+1];
T=Array[{0}&,n/Length[B]];
g=1;

For[i=1, i>0,
If[Length[U[[i]]]==n,  Tile[[g]]=T[[i]]⋃{};maxperiod[[jj,g]]=99;
For[ff=Length[Divisors[n]]-1, ff>0,
If[Mod[((Divisors[n][[ff]])+Tile[[g]]),n]⋃{}==Tile[[g]],Tile[[g]]=0;maxperiod[[jj,g]]=Divisors[n][[ff]];ff=0;
 ];
ff=ff-1];
 g=g+1;T[[i]]=T[[i-1]];mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];i=i-1,
For[j=1, j≤n,k=100;
If[Not[MemberQ[U[[i]],Mod[j,n]]], k=0;
For[m=1, m≤Length[B],
If[U[[i]]⋂ Mod[(j-B[[m]])+B,n]=={}, k=k+1; A[[i,j,m]]=Mod[j-B[[m]],n]];
m++]];
R[[j]]=k;
j++];
min=Min[R]; 
If[(Not[min==0]),
t=Catch[Do[If[R[[ii]]==min, Throw[ii]],{ii,n}]];
If[mark[[i]]≤Length[B],
For[b=mark[[i]], b<=Length[B],mark[[i]]=b+1;z=0;
If[Not[A[[i,t,b]]==200], U[[i+1]]=U[[i]]⋃Mod[A[[i,t,b]]+B,n]; T[[i+1]]=T[[i]]⋃{A[[i,t,b]]}; b=Length[B]+1;z=1;i=i+1;];
b++];
If[z==0, mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];i=i-1];
, mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];T[[i]]=T[[i-1]];i=i-1];
, mark[[i]]=1;A[[i]]=Array[200&,{n,Length[B]}];T[[i]]=T[[i-1]];i=i-1];
];
];

Tile=Drop[Tile⋃{0},1];
Length[Tile]
Tile;

For[j=1, j≤Length[Tile], P=Tile[[j]];
For[i=2, i<=Length[P], 
If[MemberQ[Tile,Mod[P-P[[i]],n]⋃{}], Tile[[j]]=0;i=Length[P]+1];
i++];
j++]; 
Tile=Drop[Tile⋃{0},1];
If[Not[Length[Tile]==0], rr=rr+1];
Tile;
jj++];
rr

For[j=1,j≤Length[maxperiod],maxperiod[[j]]=Drop[maxperiod[[j]]⋃{0},1];
j++];
maxperiod