(*{{4,8,9},{2,3,16}};*)

 
Clear[n,B,R,Tile,A,U,mark,T,g,x];

L=Expand[Cyclotomic[2,x^(9)]Cyclotomic[3,x^(8)]];

n=72;

B=Array[0&,Length[L]];
For[j=2, j≤Length[L], B[[j]]=Mod[L[[j,2]],n];
j++];
B⋃{}
(*{{4,8,9},{2,3,16}};*)

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

n=72;

B={0,8,9,16,17,25};

Length[B];

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

(*Next we shift the elements of members of Tile by 72, and evaluate a relevant entry of the Fourier transform. If it is zero, we keep the Tile, if it is non-zero we discard the Tile. The reason: if it is non-zero, then it must be zero for any tiling complement, and so the tiling complement must be periodic. *)

Dynamic[{h,p,k,j,m}]
n=144;
w2=Exp[I 2 Pi /144];

KK=Tile;

Tilea=Array[0&,200000];
h=1;
For[p=1,p≤Length[KK],
For[k=0,k≤1,
For[j=0, j≤1,
For[m=0,m≤1,
For[m2=0,m2≤1,
For[m3=0,m3≤1,
For[m4=0,m4≤1,
For[m5=0,m5≤1,
For[m6=0,m6≤1,
For[m7=0,m7≤1,
For[m8=0,m8≤1,
For[m9=0,m9≤1,
Tilea[[h]]={0,72k+KK[[p,2]],72j+KK[[p,3]],72m+KK[[p,4]],72m2+KK[[p,5]],72m3+KK[[p,6]], 72m4+KK[[p,7]],72m5+KK[[p,8]],  72m6+KK[[p,9]],72m7+KK[[p,10]], 72m8+KK[[p,11]], 72m9+KK[[p,12]]}⋃{};P=Tilea[[h]]; 
If[Not[Chop[N[1+w2^P[[2]]+w2^P[[3]]+w2^P[[4]]+w2^P[[5]]+w2^P[[6]]+w2^P[[7]]+w2^P[[8]]+w2^P[[9]]+w2^P[[10]]+w2^P[[11]]+w2^P[[12]]]]==0], Tilea[[h]]=0,h=h+1];
m9++];
m8++];
m7++];
m6++];
m5++];
m4++];
m3++];
m2++];
m++];
j++];
k++];
p++];

Tilea=Drop[Tilea⋃{0},1]
Length[Tilea]

(* It turns out that all the remaining tlies are periodic. *)

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

(* Therefore this part below is not needed, with the Fourier-selection and all that.

Dynamic[{j,v,vvv,vvvv}]

T2=Tilea;

Clear[H];

Length[T2]

n=144;
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

n=144;
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;


n=144;
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=MT;
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}];
nonperiodic=Array[0&,Length[J]];

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;
If[maxperiod[[jj,g]]==99, nonperiodic[[jj]]=J[[jj]];jj=jj+1,
If[(g>1)∧(Not[maxperiod[[jj,g]]==maxperiod[[jj,g-1]]]),nonperiodic[[jj]]=J[[jj]]; jj=jj+1]]];
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
nonperiodic=Drop[nonperiodic⋃{},1];
nonperiodic
For[j=1,j≤Length[maxperiod],maxperiod[[j]]=Drop[maxperiod[[j]]⋃{0},1];
j++];
maxperiod *)