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

(* This file is commented, the others are not (since then, i have made other files also commented). It contains most (indeed, i think all) major ideas and subroutines that are being used in any of the other files. The step-by-step method described here seems to be the most efficient that i could find. It runs very fast. However, most of the other files were written earlier and mostly contain slower algorithms because i chose non-optimal subroutines at some places. (The only optimal ones are {{3,4,9},{2,8,16}};{{3,8,9},{2,4,16}}; and this one.) 
 However, i do not want to re-write them, because i don'thave the strength anymore, so i left them in their original form. It does not mean they are not correct. They are. But sometimes a file can run for 4-5 hours, instead of a few minutes whihch would be the running time with the present optimal algorithm.*)

(* First we select the part of {{3,8,16},{2,4,9}} which is hopefully producing fewer tiles. It is usually where the LCM is larger or which has fewer elements; in this case: LCM(3,8,16)=48 (while LCM(2,4,9)=36) and both have 12 elements. 
Then apply C-M to get a universal tiling complement "B" below: *)

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

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

n=48;

B=Array[0&,Length[L]];
For[j=2, j≤Length[L], B[[j]]=Mod[L[[j,2]],n];
j++];
B⋃{}

(* The point is that B tiles with everything that is divided by each of Phi_{3,8,16}. Next we list out all complements of B mod LCM=48. In this way we find all tiles mod 48 whihch are divided by Phi_{3,8,16}. Finding these tiles mod 48 could also be ahieved by using the structure theory of Section 4 of 
C-M. That is probably quicker, and that is what we suggest in the paper. However, i only thought of it too late, and that's why i applied the present method. It gives the same result but probably much slower. Anyway, i did not have the patience to go through the structure theory of 
C-M and write yet another routine to make things run faster. However, if you want to go up to 120, 180, 200, then you should use the structure theory. 
The routine below is the "list-out-all-complements" rotuine. It is also called the "fill-out procedure" in the paper. We need to give n and B as inputs and it gives all tiling complements of B in Z_n as output.*)

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

n=48;

B={0,3,6,9}
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 put KK=Tile, and shift all elements by 0,48, or 96 in all possible ways. Thus we generate all possible tiles T in Z_144 with prescribed SA set {3,8,16}. Then we should partition the set of these tiles according to their Fourier zero-sets (which is the same as grouping them according to cyclotomic divisors). However, calculating FT is slow, tehrefore we only calculate some relevant values of it, like FT_T(1), FT_T(2), FT_T(3), FT_T(4), FT_T(6), FT_T(12), FT_T(18), FT_T(24) and so on, for divisros of 144. This will determine the cyclotomic divisors. And we don't need to calculate all of them: e.g. FT(144/3), FT(144/8), FT(144/16) certainly give zero, and also others like 144/24, 144/48 because of the C-M condition T2. Other values like 144/2,144/4,144/9 certainly give non-zero,  because of condition T1 on the tiling complements of any T. So, only a few values need to be calculated. In this particular case there is a trick: we only calculate FT_T(2)=FT_T(144/72), becasuse: if FT_T(2) is non-zero, then for any tiling complement S of T FT_S(2)=0, so \Phi_72 must divide S(x), so then S(x) is dvisible by \Phi_9,18,36,72 (18 and 36 by condition T2) so the only way to avoid periodicity of S is if \Phi_144 does not divide it. But then \Phi_144 must divide T,and therefore Phi_16,48,144 divide T, so then T is periodic. The conclusion is that the only chance for non-periodicity is if FT_T(2)=0 *)

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

KK={{0,1,2,12,13,14,24,25,26,36,37,38},{0,1,5,12,13,17,24,25,29,36,37,41},{0,1,8,12,13,20,24,25,32,36,37,44},{0,1,11,12,13,23,24,25,35,36,37,47},{0,2,4,12,14,16,24,26,28,36,38,40},{0,2,7,12,14,19,24,26,31,36,38,43},{0,2,10,12,14,22,24,26,34,36,38,46},{0,4,5,12,16,17,24,28,29,36,40,41},{0,4,8,12,16,20,24,28,32,36,40,44},{0,4,11,12,16,23,24,28,35,36,40,47},{0,5,7,12,17,19,24,29,31,36,41,43},{0,5,10,12,17,22,24,29,34,36,41,46},{0,7,8,12,19,20,24,31,32,36,43,44},{0,7,11,12,19,23,24,31,35,36,43,47},{0,8,10,12,20,22,24,32,34,36,44,46},{0,10,11,12,22,23,24,34,35,36,46,47}};

Tilea=Array[0&,200000];
h=1;
For[p=1,p≤Length[KK],
For[k=0,k≤2,
For[j=0, j≤2,
For[m=0,m≤2,
For[m2=0,m2≤2,
For[m3=0,m3≤2,
For[m4=0,m4≤2,
For[m5=0,m5≤2,
For[m6=0,m6≤2,
For[m7=0,m7≤2,
For[m8=0,m8≤2,
For[m9=0,m9≤2,
Tilea[[h]]={0,48k+KK[[p,2]],48j+KK[[p,3]],48m+KK[[p,4]],48m2+KK[[p,5]],48m3+KK[[p,6]], 48m4+KK[[p,7]],48m5+KK[[p,8]],  48m6+KK[[p,9]],48m7+KK[[p,10]], 48m8+KK[[p,11]], 48m9+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]

(* Now we apply the subroutine below to select those tiles which are non-periodic; also, it leaves alive only one representtative of rotated copies of the same tile. 
*)

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
 


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

(* In this way 18000 tiles T remain that are still possible candidates, because we found that FT_T(2)=0, and they are non-periodic. Now we partition this set of T's according to their Fourier-zero set. More precisely, we just keep one representative for each occuring zero-set. MT contains this selection*)


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

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

(* Next, we see what those Fourier-zero sets are, i.e.which cyclotomic factors they have. *)

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

(* We see that the remaining 3 zero-sets are indeed candidates for non-periodic tilings. Therefore we need to grab a representative for them (we have them at hand, above in MT) and list out all complements for each set J in MT, and see whether any complement is non-periodic. This is done in an intelligent way. First we see what the GCD of J is. For example in J=J[[2]]={0,16,24,36,50,60,74,86,100,110,124,136} below the GCD is 2. Therefore we consider J[[2]]/2 in Z_72 and list out all complements S of J[[2]]/2, and see what the maximal period of each S is (if it is 72, then S is non-periodic, and we take note of that). Clearly, if all S have the same maximal period (notequal 72), then we have no chance of producing a non-periodic tiling of Z_144 with J. If some of the maxperiods differ, we take note of that, as it gives a chance of producing non-periodic tiles with J. It does not happen in this particular case. 
Therefore we arrive at the conclusion that there is no non-periodic tiling for this partition of SA's, so we are done.*)


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

n=144;

J={{0,14,24,36,49,60,73,85,98,109,122,134},{0,16,24,36,50,60,74,86,100,110,124,136},{0,20,24,36,52,60,76,88,104,112,128,140}};
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