######################################################################
##BFF.txt: Save this file as  BFF.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read BFF.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  Nov. 12, 2015`):
print(` This is BFF.txt `):
print(`It is a package that accompanies the article `):
print(`Searching for Disjoint Covering Systems with Precisely One Repeated Modulus`):
print(`by Shalosh B. Ekhad, Aviezri  S. Fraenkel, and Doron Zeilberger `):
print(``):
print(`The package's name, BFF,  stands for  Marc Berger, Alex Felzenbaum, and Aviezri Fraenkel. `):
print(``):
print(`The article is  available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  DU, Hagdel, Halbesh, Halbesh1, Hakhlel, Hems,  `):
 print(` ImpliedDCS, IsDisj, PrimeTuples, StillAv, TovS, TovS1, Trad `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: DCS, DCS1, DCSg, DCSpc, DCSv `):
 print(` `):


elif nops([args])=1 and op(1,[args])=DCS then
print(`DCS(M,N0): inputs a positive integer M and another positive integer N0`):
print(`outputs a list of length M whose i-th item is the set of`):
print(` "almost" distinct disjoint covering systems, where all the moduli are disjoint, and distinct, except for `):
print(` the largest one, that is repeated i times. We  only list those that have the property that `):
print(` the largest modulus is the least common multiple of all the moduli, and it is <=N0. Try: `):
print(`DCS(12,100); `):

elif nops([args])=1 and op(1,[args])=DCS1 then
print(`DCS1(M,N0): inputs a positive integer M and another positive integer N0`):
print(`outputs a list of length M whose i-th item is the set of`):
print(` "almost" distinct disjoint covering systems, where all the moduli are disjoint, and distinct, except for `):
print(` the largest one, that is repeated i times. We  only list those that have the property that `):
print(` the largest modulus is the least common multiple of all the moduli, and it equals EXACTLY N0. Try: `):
print(`DCS1(12,100); `):

elif nops([args])=1 and op(1,[args])=DCSg then
print(`DCSg(M,N0,p): inputs a positive integer M and another positive integer N0`):
print(`outputs a list of length M whose i-th item is the set of`):
print(`virgin "almost" distinct disjoint covering systems, where all the moduli are disjoint, and distinct, except for`):
print(`the largest one, that is repeated i times. `):
print(`then it also outputs the implied ones, and finally all the abstract renditions of the virgin ones.`):
print(` Try: `);
print(` DCSg(12,100,p); `):

elif nops([args])=1 and op(1,[args])=DCSpc then
print(`DCSpc(): The pre-computed value of DCS(24,300);. Try `):
print(`DCSpc(); `):

elif nops([args])=1 and op(1,[args])=DCSv then
print(`DCSv(M,N0): Verbose version of DCS(M,N0) (q.v.). Try:`):
print(`DCSv(12,100);`):



elif nops([args])=1 and op(1,[args])=DU then
print(`DU(M): inputs a list of positive integers M=[m1,m2,...,mk] , and outputs the set of all lists  [a1,...,ak]`):
print(`such that the residue classes a1(mod m1), ..., ak(mod mk) are PAIRWISE disjoint. Try:`):
print(`DU([4,6]);`):

elif nops([args])=1 and op(1,[args])=Hagdel then
print(`Hagdel(Z,N): Inputs a pair Z=[a,b] where a is a divisor of N, and 0<=b<a. Outputs the subset of {0,...,N-1}`):
print(`occupied by the resiude class b (mod a). Try:`):
print(`Hagdel([7,4],28); `):


elif nops([args])=1 and op(1,[args])=Hakhlel then
print(`Hakhlel(L,p): Given a concrete DCS [L1,L2, ..., [m,timesrepeated]]`):
print(`finds the infinite family implied by it, phrased in terms of p[1],p[2], ..., followed by the primes participating.`):
print(` Try:Hakhlel([3, 15, [30, 18]],p); `):

elif nops([args])=1 and op(1,[args])=Halbesh then
print(`Halbesh(L,N): Inputs a list L of divisors of L outputs a list of pairs [a,b] that`):
print(`the residue classes b(mod a) do not intersect, the last item N with the`):
print(`remaining residue classes. `):
print(`Halbesh([3,6],12);`):

elif nops([args])=1 and op(1,[args])=Halbesh1 then
print(`Halbesh1(L,N): Inputs a list L of divisors of L outputs a list of pairs [a,b] that`):
print(`the residue classes b(mod a) do not intersect. It is picked randomly. Try:`):
print(`Halbes1([3,6],12);`):

elif nops([args])=1 and op(1,[args])=Hems then
print(`Hems(L,N): all the legal continuations of the sequence L of divisors of N. Try:`):
print(`Hems([3],18);`):

elif nops([args])=1 and op(1,[args])=ImpliedDCS then
print(`ImpliedDCS(A,Ma,p): given an abstract DCS A phrased in terms of p, outputs all the implied concrete `):
print(`almost distinct DCSs`):
print(`that have the multiplicity of the largest part<=Ma. Try:`):
 print(` ImpliedDCS([[p[1],[p[1]^2,p[1]^2-p[1]]],[2]]  ,30,p); `):

elif nops([args])=1 and op(1,[args])=IsDisj then
print(`IsDisj(P1,P2): inputs two pairs P1=[m1,a1], P2=[m2,a2] returns true iff (a1 mod m1) and (a2 mod m2) are disjoint. Try:`):
print(`IsDisj([3,2],[15,3]);`):


elif nops([args])=1 and op(1,[args])=PrimeTuples then
print(`PrimeTuples(P,M): all the prime tuples of length nops(P) that are bigger than P and less than M.`):
print(` Try: `):
print(` PrimeTuples([2,3,5],[11,11,11]); `):

elif nops([args])=1 and op(1,[args])=StillAv then
print(`StillAv(L,N,i): Inputs a list of L, of pairs [a,b], where a is a divisor of N and b is a residue class mod a (0<=b<a)`):
print(`and i is another divisor of N, outputs the set of residue classes that i can have without conflicting with the`):
print(`already established residue classes. Try:`):
print(`StillAv([[3,1],[6,0]],36,9);`):


elif nops([args])=1 and op(1,[args])=TovS then
print(`TovS(N): the set of sequences of divisors of N of any length  that have density less than 1 and with the`):
print(`property that gcd(i,j)>1 for each of them. `):
print(`the output is a pair. The first output is the set of m's,and the second is a table such that`):
print(`T[m] is the set of good sequences of divisors of N. Try:`):
print(`TovS(36); `):



elif nops([args])=1 and op(1,[args])=TovS1 then
print(`TovS1(N,k): the set of sequences of divisors of N of length k that have density less than 1 and with the`):
print(`property that gcd(i,j)>1 for each of them. Try:`):
print(`TovS1(36,3); `):
print(``):

elif nops([args])=1 and op(1,[args])=Trad then
print(`Trad(L,P,p): given a list L [[p1,m1],[p2,m2], ...] of concrete primes P[1],P[2] with multiplicities,`):
print(`outputs the expression in terms of abstract primes,`):
print(` Try: `):
print(` Trad([[3,1],[5,3]],[3,5],p); `):


else
print(`There is no ezra for`,args):
fi:
 
end:


with(numtheory):


#Hems(L,N): all the legal continuations of the sequence L of divisors of N. Try
#Hems([3],24);
Hems:=proc(L,N) local i1,lu,gu,gadol,ku,lu1:

if L=[] then
 RETURN(FAIL):
fi:

ku:=add(1/L[i1],i1=1..nops(L)):

if ku>=1 then
   RETURN({}):
fi:

lu:=divisors(N) minus {1,2,N}:

gadol:=L[nops(L)]:

gu:={}:

for lu1 in lu do
 if lu1>gadol and min(seq(gcd(L[i1],lu1),i1=1..nops(L)))>1 then
   if ku+1/lu1<1 then
     gu:=gu union {[op(L),lu1]}:
   fi:
 fi:
od:

gu:

end:


#TovS1(N,k): the set of sequences of divisors of N of length k that have density less than 1 and with the
#property that gcd(i,j)>1 for each of them

TovS1:=proc(N,k) local gu,gu1, lu,i1:
option remember:

lu:=divisors(N) minus {1,2,N}:

if k=1 then
 RETURN({seq([lu[i1]],i1=1..nops(lu))}):
fi:

gu:=TovS1(N,k-1):

{seq(op(Hems(gu1,N)),gu1 in gu)}:

end:

#Hagdel(Z,N): Inputs a pair Z=[a,b] where a is a divisor of N, and 0<=b<a. Outputs the subset of {0,...,N-1}
#occupied by the resiude class b (mod a). Try:
#Hagdel([7,4],28);
Hagdel:=proc(Z,N) local a,b,j:

if not (type(Z,list) and nops(Z)=2 and type(Z[1],integer) and  type(Z[2],integer) and type(N, integer)
and type(N/Z[1],integer) and Z[1]>0 and Z[2]>=0 and Z[2]<Z[1]) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

a:=Z[1]: b:=Z[2]:
{seq(b+j*a,j=0..N/a-1)}:

end:

#StillAv(L,N,i): Inputs a list of L, of pairs [a,b], where a is a divisor of N and b is a residue class mod a (0<=b<a)
#and i is another divisor of N, outputs the set of residue classes that i can have without conflicting with the
#already established residue classes. Try:
#StillAv([[3,1],[6,0]],36,9);
StillAv:=proc(L,N,i) local gu,i1,tafus:

if not (type(L,list) and type(N,integer) and N>0 and type(N/i,integer) ) then
 RETURN(FAIL):
fi:

tafus:={seq(op(Hagdel(L[i1],N)),i1=1..nops(L))}:


gu:={}:

for i1 from 0 to i-1 do


 if Hagdel([i,i1],N) intersect tafus={} then
  gu:=gu union {i1}:
 fi:
od:
gu:


end:


#Halbesh1Old(L,N): Inputs a list L of divisors of L outputs a list of pairs [a,b] that
#the residue classes b(mod a) do not intersect. It is picked randomly. Try
#Halbes1Old([3,6],12);
Halbesh1Old:=proc(L,N) local i,L1,gu,kha:

if nops(L)=0 then
 RETURN(FAIL):
fi:

L1:=[[L[1],0]]:

for i from 2 to nops(L) do
gu:=StillAv(L1,N,L[i]) :
if gu={} then
 RETURN(FAIL):
else
kha:=gu[1]:
L1:=[op(L1),[L[i],kha]]:
fi:
od:

L1:

end:




#TovS(N): the set of sequences of distinct legal divisors of N such that the density is less than 1.
#the output is a pair. The first output is the set of m's,and the second is a table such that
#T[m] is the set of good sequences.
#Try:
#TovS(12);
#TOvS(36);
TovS:=proc(N) local gu,gu1, k,m,lu,lu1,i,i1,ru,T,ru1:

gu1:=TovS1(N,1):

gu:=gu1:

for k from 1 while gu1<>{} do
gu1:=TovS1(N,k):
gu:=gu union gu1:
od:

lu:={}:
ru:={}:
for i from 1 to nops(gu) do
m:=N-add(N/gu[i][i1],i1=1..nops(gu[i]) ):
lu:=lu union {[gu[i],m]}:
ru:=ru union {m}:
od:

for ru1 in ru do
 T[ru1]:={}:
od:

for lu1 in lu do
 T[lu1[2]]:=T[lu1[2]] union {lu1[1]}:
od:

ru,op(T):

end:


#Halbesh(L,N): Inputs a list L of divisors of L outputs a list of pairs [a,b] that
#the residue classes b(mod a) do not intersect, the last item N with the
#remaining residue classes. 
#Halbes([3,6],12);
Halbesh:=proc(L,N) local gu,mu:
gu:=Halbesh1(L,N):

if gu=FAIL then
  RETURN(FAIL):
fi:

mu:=StillAv(gu,N,N):

[op(gu),[N,mu]]:

end:



DCSpc:=proc():
[{}, {}, {}, {[3, [6, 4]]}, {}, {[3, [9, 6]], [4, [8, 6]], [3, 6, [12, 6]]}, {[
4, 6, [12, 7]], [3, 6, 9, [18, 7]]}, {[3, [12, 8]], [5, [10, 8]]}, {[4, [12, 9]
], [3, 6, [18, 9]], [4, 6, 8, 12, [24, 9]], [3, 6, 9, 12, 18, [36, 9]]}, {[3, [
15, 10]], [6, [12, 10]], [3, 9, [18, 10]], [4, 8, [16, 10]], [3, 6, 12, [24, 10
]]}, {[4, 6, 8, [24, 11]], [3, 6, 9, 12, [36, 11]]}, {[3, [18, 12]], [4, [16, 
12]], [5, [15, 12]], [7, [14, 12]], [3, 6, [24, 12]], [4, 6, 12, [24, 12]], [3,
6, 9, 18, [36, 12]]}, {[4, 10, [20, 13]], [6, 9, [18, 13]], [3, 6, 15, [30, 13]
], [4, 8, 12, [24, 13]], [3, 6, 12, 18, [36, 13]], [4, 6, 8, 12, 16, 24, [48, 
13]], [3, 6, 9, 12, 18, 24, 36, [72, 13]]}, {[3, [21, 14]], [8, [16, 14]], [3,
12, [24, 14]], [4, 6, [24, 14]], [5, 10, [20, 14]], [3, 6, 9, [36, 14]]}, {[4,
[20, 15]], [6, [18, 15]], [3, 6, [30, 15]], [3, 9, [27, 15]], [4, 8, [24, 15]],
[3, 6, 12, [36, 15]], [6, 8, 12, [24, 15]], [3, 9, 12, 18, [36, 15]], [4, 6, 8,
12, 16, [48, 15]], [3, 6, 9, 12, 18, 24, [72, 15]]}, {[3, [24, 16]], [5, [20, 
16]], [9, [18, 16]], [4, 12, [24, 16]], [3, 6, 18, [36, 16]], [4, 6, 12, 18, [
36, 16]], [3, 6, 9, 18, 27, [54, 16]], [4, 6, 8, 12, 24, [48, 16]], [3, 6, 9, 
12, 18, 36, [72, 16]]}, {[6, 8, [24, 17]], [3, 9, 12, [36, 17]], [4, 6, 8, 16,
24, [48, 17]], [3, 6, 9, 12, 24, 36, [72, 17]]}, {[3, [27, 18]], [4, [24, 18]],
[7, [21, 18]], [10, [20, 18]], [3, 6, [36, 18]], [3, 15, [30, 18]], [6, 12, [24
, 18]], [3, 9, 18, [36, 18]], [4, 6, 12, [36, 18]], [4, 8, 16, [32, 18]], [3, 6
, 9, 18, [54, 18]], [3, 6, 12, 24, [48, 18]], [4, 6, 8, 12, [48, 18]], [3, 6, 9
, 12, 18, [72, 18]], [4, 6, 8, 12, 18, 24, 36, [72, 18]], [3, 6, 9, 12, 18, 27,
36, 54, [108, 18]]}, {[4, 14, [28, 19]], [8, 12, [24, 19]], [3, 6, 21, [42, 19]
], [3, 12, 18, [36, 19]], [4, 6, 18, [36, 19]], [5, 10, 15, [30, 19]], [3, 6, 9
, 27, [54, 19]], [4, 6, 8, 16, [48, 19]], [4, 8, 10, 20, [40, 19]], [3, 6, 9, 
12, 24, [72, 19]], [3, 6, 12, 15, 30, [60, 19]], [4, 6, 12, 16, 24, [48, 19]],
[3, 6, 9, 18, 24, 36, [72, 19]]}, {[3, [30, 20]], [5, [25, 20]], [6, [24, 20]],
[11, [22, 20]], [3, 9, [36, 20]], [4, 8, [32, 20]], [3, 6, 12, [48, 20]], [6, 
10, 15, [30, 20]], [4, 6, 8, 24, [48, 20]], [3, 6, 9, 12, 36, [72, 20]], [4, 6,
8, 12, 18, 24, [72, 20]], [3, 6, 9, 12, 18, 27, 36, [108, 20]]}, {[4, [28, 21]]
, [8, [24, 21]], [3, 6, [42, 21]], [3, 12, [36, 21]], [4, 6, [36, 21]], [5, 10,
[30, 21]], [3, 6, 9, [54, 21]], [4, 8, 10, [40, 21]], [3, 6, 12, 15, [60, 21]],
[4, 6, 12, 16, [48, 21]], [6, 9, 12, 18, [36, 21]], [3, 6, 9, 18, 24, [72, 21]]
, [4, 8, 12, 16, 24, [48, 21]], [3, 6, 12, 18, 24, 36, [72, 21]], [4, 6, 8, 12,
18, 36, [72, 21]], [3, 6, 9, 12, 18, 27, 54, [108, 21]], [4, 6, 8, 12, 16, 24,
32, 48, [96, 21]], [3, 6, 9, 12, 18, 24, 36, 48, 72, [144, 21]]}, {[3, [33, 22]
], [12, [24, 22]], [3, 18, [36, 22]], [4, 16, [32, 22]], [5, 15, [30, 22]], [6,
10, [30, 22]], [7, 14, [28, 22]], [3, 6, 24, [48, 22]], [3, 9, 15, [45, 22]], [
4, 6, 8, [48, 22]], [4, 12, 18, [36, 22]], [3, 6, 9, 12, [72, 22]], [3, 6, 18,
27, [54, 22]], [4, 6, 12, 24, [48, 22]], [3, 6, 9, 18, 36, [72, 22]], [4, 6, 8,
12, 24, 36, [72, 22]], [3, 6, 9, 12, 18, 36, 54, [108, 22]], [4, 6, 8, 12, 16,
18, 24, 36, 48, 72, [144, 22]], [3, 6, 9, 12, 18, 24, 27, 36, 54, 72, 108, [216
, 22]]}, {[6, 15, [30, 23]], [4, 8, 20, [40, 23]], [6, 9, 12, [36, 23]], [3, 6,
12, 30, [60, 23]], [4, 6, 16, 24, [48, 23]], [4, 8, 12, 16, [48, 23]], [3, 6, 9
, 24, 36, [72, 23]], [3, 6, 12, 18, 24, [72, 23]], [4, 6, 8, 12, 18, [72, 23]],
[3, 6, 9, 12, 18, 27, [108, 23]], [4, 6, 8, 12, 16, 24, 32, [96, 23]], [3, 6, 9
, 12, 18, 24, 36, 48, [144, 23]]}, {[3, [36, 24]], [4, [32, 24]], [5, [30, 24]]
, [7, [28, 24]], [9, [27, 24]], [13, [26, 24]], [3, 6, [48, 24]], [4, 12, [36,
24]], [3, 6, 18, [54, 24]], [4, 6, 12, [48, 24]], [4, 10, 20, [40, 24]], [6, 9,
18, [36, 24]], [3, 6, 9, 18, [72, 24]], [3, 6, 15, 30, [60, 24]], [4, 8, 12, 24
, [48, 24]], [3, 6, 12, 18, 36, [72, 24]], [4, 6, 8, 12, 24, [72, 24]], [3, 6,
9, 12, 18, 36, [108, 24]], [4, 6, 8, 18, 24, 36, [72, 24]], [3, 6, 9, 12, 27, 
36, 54, [108, 24]], [4, 6, 8, 12, 16, 24, 48, [96, 24]], [3, 6, 9, 12, 18, 24,
36, 72, [144, 24]], [4, 6, 8, 12, 16, 18, 24, 36, 48, [144, 24]], [3, 6, 9, 12,
18, 24, 27, 36, 54, 72, [216, 24]]}]:
end:

#DCS(M,N0): inputs a positive integer M and another positive integer N0
#outputs a list of length M whose i-th item is the set of
#"almost" distinct disjoint covering systems, where all the moduli are disjoint, and distinct, except for
#the largest one, that is repeated i times. We  only list those that have the property that
#the largest modulus is the least common multiple of all the moduli, and it is <=N0. Try:
#DCS(12,100);

DCS:=proc(M,N0) local m,T1,n,gu,S,kv,kv1,i, muams,muam, muam1,i1:
option remember:
for m from 1 to M do
 T1[m]:={}:
od:


for n from 1 to N0 do
 gu:=TovS(n):
 kv:=gu[1] intersect {seq(i1,i1=1..M)}:
 S:=gu[2]:

  for kv1 in kv do
    muams:=S[kv1]:
     for muam in muams do
       muam1:=Halbesh(muam,n) :
       if muam1<>FAIL then
        T1[kv1]:=T1[kv1] union {[op(muam),[n,kv1]]}:
         else
                   if DU(muam)<>{} then
                   T1[kv1]:=T1[kv1] union {[op(muam),[n,kv1]]}:
               fi:
       fi:
   od:
  od:
od:

[seq(T1[i],i=1..M)]:


end:


#DCS1(M,N0): inputs a positive integer M and another positive integer N0
#outputs a list of length M whose i-th item is the set of
#"almost" distinct disjoint covering systems, where all the moduli are disjoint, and distinct, except for
#the largest one, that is repeated i times. We  only list those that have the property that
#the largest modulus is the least common multiple of all the moduli, and it equals N0. Try:
#DCS1(12,100);

DCS1:=proc(M,N0) local m,T1,n,gu,S,kv,kv1,i, muams,muam, muam1,i1:

for m from 1 to M do
 T1[m]:={}:
od:


for n from N0 to N0 do
 gu:=TovS(n):
 kv:=gu[1] intersect {seq(i1,i1=1..M)}:
 S:=gu[2]:

  for kv1 in kv do
    muams:=S[kv1]:
     for muam in muams do
       muam1:=Halbesh(muam,n) :
       if muam1<>FAIL then
        T1[kv1]:=T1[kv1] union {[op(muam),[n,kv1]]}:
         else
                   if DU(muam)<>{} then
                   T1[kv1]:=T1[kv1] union {[op(muam),[n,kv1]]}:
               fi:
       fi:
   od:
  od:
od:

[seq(T1[i],i=1..M)]:


end:


#IsDisj(P1,P2): inputs two pairs P1=[m1,a1], P2=[m2,a2] returns true iff (a1 mod m1) and (a2 mod m2) are disjoint. Try:
#IsDisj([3,2],[15,3]);
IsDisj:=proc(P1,P2) local a1,a2,m1,m2,M,i1:
m1:=P1[1]:m2:=P2[1]:a1:=P1[2]:a2:=P2[2]:
M:=lcm(m1,m2):
evalb({seq(a1+m1*i1,i1=0..M/m1)} intersect {seq(a2+m2*i1,i1=0..M/m2)}={}):
end:


#DU(M): inputs a list of positive integers M=[m1,m2,...,mk] , and outputs the set of all lists  [a1,...,ak]
#such that the residue classes a1(mod m1), ..., ak(mod mk) are PAIRWISE disjoint. Try:
#DU([4,6]);
DU:=proc(L) local gu,i1,L1,mu,mu1,i2,n,cu:
option remember:

n:=nops(L):

if n=0 then
 RETURN(FAIL):
fi:

if n=1 then
RETURN({seq([[L[1],i1]],i1=0..L[1]-1)}):
fi:

L1:=[op(1..n-1,L)]:

mu:=DU(L1):
cu:=L[n]:

gu:={}:

for mu1 in mu do

 for i1 from 0 to cu-1 do
  
  if {seq(IsDisj(mu1[i2],[cu,i1]) ,i2=1..nops(mu1))}={true} then
    gu:=gu union {[op(mu1),[cu,i1] ] }:
  fi:
 od:
od:
gu:

end:



#DCSv(M,N0): Verbose version of DCS(M,N0) (q.v.). Try:
#DCSv(12,100);

DCSv:=proc(M,N0) local t0,gu,mu,j,mu1,i:

t0:=time():

gu:=DCS(M,N0):

print(``):
print(`All the Exact Covering Systems with Distinct Moduli Except for the Largest that is repeated Up to`, M,  ` Times `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`An  EXACT covering System is a set of congruences a_i(mod m_i), i=1..k, with m_1<=...<=m_k such that every positive integer`):
print(` belongs to exactly one of these residue classes. A classical theorem of Erdos, Davenport, Newman, and Mirsky`):
print(`asserts that the largest modulus must be repeated at least twice.`):
print(``):
print(`We are interested in such systems where all the moduli are distinct, except the largest, that is repeated a`):
print(`specified number of times. `):
print(`In other words, we have m_1<m_2<...<m_{k-1} followed by m_k repeated. `):
print(``):
print(`In this exciting article, we will find ALL such systems where the smallest modulus is at least 3, and the largest modulus`):
print(` is repeated exactly M times, for M from 2 to`, M, `    . ` ):
print(`We only list those systems where the largest modulus is the least common multiple of all the moduli, and it is`):
print(`at most`, N0):

print(``):
print(`The complete list up to 9 repetitions was done in:`):
print(` M.A. Berger, A. Felzenbaum, and A. Fraenkel, Acta Arith. (1988), 171-182`):
print(`This was extended to up to 12 repetitions by  M. Zeleke and J. Simpson`):
print(`in Adv. Appl. Math. 23(1999), 322-332. `):
print(``):

print(`We  exclude the utterly trivial systems that only have one modulus. `):
print(``):
print(` Since `, N0 , `is large enough, we conjecture that this covers all possibilities.`):
print(``):

for i from 1 to M do
print(``):

print(`--------------------------------------------------------------------`):
mu:=gu[i]:
if nops(mu)=0 then
  print(`There are no systems where the largest modulus is repeated  `, i, `times . `):

elif nops(mu)=1 then
print(`There is exactly one system where the largest modulus is repeated  `, i, `times . Here it is `):
print(``):
mu1:=mu[1]:
print(`The moduli are:`, op(1..nops(mu1)-1, mu1), `followed by`, mu1[nops(mu1)][1], `that is repeated`, i, ` times . `):
else
print(`There are exactly`, nops(mu), `systems where the largest modulus is repeated  `, i, `times . Here they are.`):
print(``):

for j from 1 to nops(mu) do
mu1:=mu[j]:
print(`The moduli are:`, op(1..nops(mu1)-1, mu1), `followed by`, mu1[nops(mu1)][1],`that is repeated`, i, ` times . `):
od:
fi:
od:
print(``):
print(`----------------------------------------------------------------`):
print(``):

print(`This ends this article, that took`, time() - t0, `seconds. to generate. `):

end:



#Halbesh1(L,N): Inputs a list L of divisors of L outputs a list of pairs [a,b] that
#the residue classes b(mod a) do not intersect. It is picked randomly. Try
#Halbes1([3,6],12);
Halbesh1:=proc(L,N) local gu,i:

for i from 1 to 20 do
 gu:=Halbesh11(L,N):
  if gu<>FAIL then
    RETURN(gu):
  fi:
od:
FAIL:
end:


#Halbesh11(L,N): Inputs a list L of divisors of L outputs a list of pairs [a,b] that
#the residue classes b(mod a) do not intersect. It is picked randomly. Try
#Halbes11([3,6],12);
Halbesh11:=proc(L,N) local i,L1,gu,kha:

if nops(L)=0 then
 RETURN(FAIL):
fi:

L1:=[[L[1],0]]:

for i from 2 to nops(L) do
gu:=StillAv(L1,N,L[i]) :
if gu={} then
 RETURN(FAIL):
else
kha:=gu[rand(1..nops(gu))()]:
L1:=[op(L1),[L[i],kha]]:
fi:
od:

L1:

end:




#Trad(L,P,p): given a list L [[p1,m1],[p2,m2], ...] of concrete primes P[1],P[2] with multiplicities,
#outputs the expression in terms of abstract primes,
#Try:
#Trad([[3,1],[5,3]],[3,5],p);
Trad:=proc(L,P,p) local k,i,P1,M1:
k:=nops(P):
P1:=[seq(L[i][1],i=1..nops(L))]:
M1:=[seq(L[i][2],i=1..nops(L))]:
P1:=subs({seq(P[i]=p[i],i=1..k)},P1):
mul(P1[i]^M1[i],i=1..nops(L)):
end:


#Hakhlel(L,p): Given a concrete DCS [L1,L2, ..., [m,timesrepeated]], and a symbol p,
#finds the infinite family implied by it in terms of abstract primes p[1],p[2], ..., followed by the primes participating.
#Try:Hakhlel([3, 15, [30, 18]],p);
Hakhlel:=proc(L,p) local L1,N,gu,P,i,M:
N:=L[nops(L)][1]:
gu:=ifactors(N)[2]:
P:=[seq(gu[i][1],i=1..nops(gu))];
L1:=[op(1..nops(L)-1,L),N]:
L1:=[seq(ifactors(L1[i])[2],i=1..nops(L1))]:
L1:=[seq(Trad(L1[i],P,p),i=1..nops(L1))]:
N:=L1[nops(L1)]:
M:=normal((1-add(1/L1[i],i=1..nops(L1)-1))*N):

[[op(1..nops(L1)-1,L1),[N,M]],P]:
end:

#PrimeTuples(P,M): all the prime tuples of length nops(P) that are bigger than P and less than M.
#Try:
#PrimeTuples([2,3,5],11);
PrimeTuples:=proc(P,M) local P1,M1,ru,p,gu,gu1,ru1,vu,vu1:
option remember:

if P=[] then
 RETURN({[]}):
fi:

P1:=[op(1..nops(P)-1,P)]:
M1:=[op(1..nops(M)-1,M)]:

gu:=PrimeTuples(P1,M1):

ru:={}:

p:=P[nops(P)]:

while p<=M[nops(M)] do
 ru:=ru union {p}:
 p:=nextprime(p):
od:

vu:={seq(seq([op(gu1),ru1],gu1 in gu), ru1 in ru)}:

gu:={}:

for vu1 in vu do
 if nops(convert(vu1,set))=nops(vu1) then
  gu:=gu union {vu1}:
 fi:
od:
gu:
end:

#ImpliedDCS(A,Ma,p): given an abstract DCS A phrased in terms of p, outputs all the implied concrete 
#almost distinct DCSs
#that have the multiplicity of the largest part<=Ma.
#Try: ImpliedDCS([[p[1],[p[1]^2,p[1]^2-p[1]]],[2]]  ,30,p);
ImpliedDCS:=proc(A,Ma,p) local k,x,A1,P,i,x1,lu,lu1,KAMA,M,gu:
A1:=A[1]:
P:=A[2]:
k:=nops(P):
KAMA:=A1[nops(A1)][2]:



gu:=subs({seq(p[i]=x,i=1..k)},KAMA):



for x1 from 1 while subs(x=x1,gu)<=Ma+30 do od:
 M:=x1+1:


gu:={}:

lu:=PrimeTuples(P,[M$k]):

for lu1 in lu do

 if subs({seq(p[i]=lu1[i],i=1..k)},KAMA)<=Ma then
  gu:=gu union  {subs({seq(p[i]=lu1[i],i=1..k)},A1)}:
fi:
od:

gu:


end:


#DCSg(M,N0,p): inputs a positive integer M and another positive integer N0
#outputs a list of length M whose i-th item is the set of
#virgin "almost" distinct disjoint covering systems, where all the moduli are disjoint, and distinct, except for
#the largest one, that is repeated i times. 
#then it also outputs the implied ones, and finally all the abstract renditions of the virgin ones.
#Try:
#DCSg(12,100,p);

DCSg:=proc(M,N0,p) local gu,IM,VIR,i,IMP,ABS,ku, kv,kv1,ru,ru1:
option remember:
gu:=DCS(M,N0):

IMP:={}:

for i from 1 to M do
 IM[i]:={}:
 VIR[i]:={}:
 ABS[i]:={}:
od:




for i from 1 to M do
 kv:=gu[i]:

 for kv1 in kv do
   if not member(kv1,IMP) then
    VIR[i]:=VIR[i] union {kv1}:
    ku:=Hakhlel(kv1,p):
    ABS[i]:= ABS[i] union {ku}:
    ru:=ImpliedDCS(ku,M,p) minus {kv1}:

    IMP:=IMP union ru :
     for ru1 in ru do
      IM[ru1[nops(ru1)][2]]:=     IM[ru1[nops(ru1)][2]] union {ru1}:
     od:   fi:
 od:
od:  



[seq(VIR[i],i=1..M)],  [seq(IM[i],i=1..M)], [seq(ABS[i],i=1..M)]:

end:



#DCSgv(M,N0,p): Verbose version of DCSg(M,N0,p) (q.v.). Try:
#DCSgv(12,100,p);

DCSgv:=proc(M,N0,p) local t0,gu,mu1, mu2, j,i,k,i1,mu11,mu31:

t0:=time():

gu:=DCSg(M,N0,p):

print(``):
print(`All the Exact Covering Systems with Distinct Moduli Except for the Largest that is repeated Up to`, M,  ` Times `):
print(`And the corresponding abstract families `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`An  EXACT covering System is a set of congruences a_i(mod m_i), i=1..k, with m_1<=...<=m_k such that every positive integer`):
print(` belongs to exactly one of these residue classes. A classical theorem of Erdos, Davenport, Newman, and Mirsky`):
print(`asserts that the largest modulus must be repeated at least twice.`):
print(``):
print(`We are interested in such systems where all the moduli are distinct, except the largest, that is repeated a`):
print(`specified number of times. `):
print(`In other words, we have m_1<m_2<...<m_{k-1} followed by m_k repeated. `):
print(``):
print(`In this exciting article, we will find ALL such systems where the smallest modulus is at least 3, and the largest modulus`):
print(` is repeated exactly M times, for M from 2 to`, M, `    . ` ):

print(``):
print(`The complete list up to 9 repetitions was done in:`):
print(` M.A. Berger, A. Felzenbaum, and A. Fraenkel, Acta Arith. (1988), 171-182`):
print(`This was extended to up to 12 repetitions by  M. Zeleke and J. Simpson`):
print(`in Adv. Appl. Math. 23(1999), 322-332. `):
print(``):

print(`We  exclude the utterly trivial systems that only have one modulus. `):
print(`We only list those systems where the largest modulus is the least common multiple of all the moduli, and it is`):
print(`at most`, N0):
print(``):
print(` Since `, N0 , `is large enough, we conjecture that this covers all possibilities.`):
print(``):

print(``):
print(``):
print(`For each number of repetitions, we list the primitive ones, that start an infinite `):
print(`family, and the corresponding abstract family that they generate, and also the ones`):
print(`that are implied from an already existing family`):

for i from 1 to M do
print(``):

print(`--------------------------------------------------------------------`):

mu1:=gu[1][i]:mu2:=gu[2][i]: 

if nops(mu1)=0 and nops(mu2)=0 then
  print(`There are no systems where the largest modulus is repeated  `, i, `times . `):
fi:

if mu1<>{} then
print(`There are`, nops(mu1), `primitive systems where the moduli are distinct`):
print(`except for the largest one that is repeated`, i, `times. `):
print(`Here there are, together with the infinite family of such systems`):
print(`implied by each. `):
print(``):

for j from 1 to nops(mu1) do
mu11:=mu1[j]:
print(``):
print(` The `, j, ` -th such system is `):
print(``):
print(`The moduli are:`, op(1..nops(mu11)-1, mu11), `followed by`, mu11[nops(mu11)][1],`that is repeated`, i, ` times . `):
print(``):

print(`This primitive system implies the infinite family`):

mu31:=Hakhlel(mu11,p):
k:=nops(mu31[2]):


print(`The moduli are:`, op(1..nops(mu31[1])-1, mu31[1]), `followed by`, mu31[1][nops(mu31[1])][1],`that is repeated`):
print( mu31[1][nops(mu31[1])][2] , `times `):
print(``):
print(` where` ,seq(p[i1],i1=1..k)):
print(`are any distinct primes satisfying`):
print(seq(p[i1]>=mu31[2][i1],i1=1..k)):
od:
fi:

if mu2<>{} then
print(`In addition there are`, nops(mu2), `systems that are implied by previously found systems. Here they are.`):


for j from 1 to nops(mu2) do
mu11:=mu2[j]:
print(``):
print(` The `, j, ` -th such system is `):
print(``):
print(`The moduli are:`, op(1..nops(mu11)-1, mu11), `followed by`, mu11[nops(mu11)][1],`that is repeated`, i, ` times . `):
print(``):

od:

fi:

od:
print(``):
print(`----------------------------------------------------------------`):
print(``):

print(`This ends this article, that took`, time() - t0, `seconds. to generate. `):

end: