################################
#This is a Maple package written 
#to help efficiently search for
#partition identities

#by Mingjia Yang

#Aug. 26, 2020
################################

print(`Search:`):
print(`A Maple package written`):
print(`to help efficiently search for partition identities`):
print(`Requires Frank Garvan's q-series package, available at`):
print(`http://www.qseries.org/fgarvan/qmaple/qseries/index.html`):
print(`For help, type Help();`):



with(combinat):
with(qseries):

Help:=proc()

if args=NULL then
 print(`The main procedures are: GP, GxnSeq,Search `):
 print(` `):

elif nops([args])=1 and op(1,[args])=GxnSeq then
print(`GxnSeq(N,A,Mod,B,IC): the first N terms of the sequence "number of `):
print(` partitions of n avoiding the partition patterns A globally,`):
print(`patterns in Mod locally (according to mod conditions),and patterns in B at the beginning,`):
print(` in addition to having no sub-partitions in IC (initial conditions)". Try: `):
print(`GxnSeq(50,{[0],[1]},[],{},[]);`):
print(`GxnSeq(50,{[0],[1]},[],{},[[1]]);`):
print(`GxnSeq(50,{[0],[1]},[{[2]},{[3]},{[2],[3]}],{},[[2]])`):

elif nops([args])=1 and op(1,[args])=Search then
print(`Search(N,A,Mod,B,IC,S): Searching for partition identities such that `):
print(`the sum side is the generating function for partitions that avoid patterns in a globally,`):
print(`[m1,m2..] according to mod conditions, ini as sub-partitions, B at the beginning`):
print(`where a is in A, ini is in powerset(IC),[m1,m2..] is in `):
print(` the cartesian product of the power sets of the members of Mod`):
print(` and prodmake outputs members in S`):
print(` Try: Search(30,{{[0],[1]}},[{[2]},{[3]},{[2],[3]}],{},[[1]],{1,0,-1,-2});`):

print(``):

elif nops([args])=1 and op(1,[args])=GP then
print(`GP(m,n,A,Mod,B,IC): the number of `):
print(` partitions of n with the largest part m, avoiding the partition patterns A globally,`):
print(`patterns in Mod locally (according to mod conditions),and patterns in B at the beginning,`):
print(` in addition to having no sub-partitions in IC (initial conditions)". Try: `):
print(`GP(8,50,{[0],[1]},[{[2]},{[3]},{[2],[3]}],{},[[2]])`):

elif nops([args])=1 and op(1,[args])=Yeled then
print(`Yeled(d,S): The child of the set of patterns S when the difference between the largest and second-largest part is d.`):
print(`Try: `):
print(`Yeled(2,{[2,1]}); `):

elif nops([args])=1 and op(1,[args])=Yeled1 then
print(`Yeled1(d,p): The child of the pattern p when the difference between the largest and second-largest part is d.`):
print(`Try: `):
print(`Yeled1(2,[2,1]); `):

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



#######################################################################
#Beginning of borrowing from RPpos.txt, available at:                 #
#https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/rpr.html#
#######################################################################

#IsKickable(p,S): given a list of prefix patterns S and p a member, can p be kicked out?
IsKickable:=proc(p,S) local p1:

for p1 in S minus {p} do
 if nops(p1)<=nops(p) and [op(1..nops(p1),p)]=p1 then
  RETURN(true):
 fi:
od:

false:

end:

#CleanUp1(S): One step  of cleaning up
CleanUp1:=proc(S) local p,S1:
S1:=S:

for p in S do
 if IsKickable(p,S) then
     S1:=S1 minus {p}:
 fi:
od:
S1:
end:

CleanUp:=proc(S) local S1,S2,S3:
S1:=S:
S2:=CleanUp1(S1):

while S1<>S2 do

S3:=CleanUp1(S2):

S2:=S3:
S1:=S2:
od:

S2:

end:

#Yeled1(d,p): The child of the pattern p when the difference between the largest and second-largest part is d.
#Try:
#Yeled1(2,[2,1]);
Yeled1:=proc(d,p)
option remember:

if d<0 or p=[] then
  RETURN(FAIL):
fi:

if p[1]=d then
 RETURN({[op(2..nops(p),p)]}):
else
 RETURN({}):
fi:

end:

#Yeled(d,S): the set of children of the set of patterns S when the difference between the largest and second-largest terms is d. Try:
#Yeled(1,{[0],[1,2],[1,3],[2]});
Yeled:=proc(d,S) local p,gu,mu:
option remember:
gu:={}:

for p in S do
  mu:=Yeled1(d,p):

  if mu<>FAIL then
    gu:=gu union mu:
  fi:

od:
CleanUp(gu):
end:

####################################
#Ending of borrowing from RPpos.txt#
####################################

                                                        
#GP(m,n,A,Mod,B,IC): the number of 
# partitions of n with largest part=m avoiding the partition patterns A globally, patterns in Mod locally (according to mod conditions),
#and patterns in B at the beginning, in addition to having no sub-partitions in IC (initial conditions)

GP:=proc(m,n,A,Mod,B,IC) local gu,m1,k,i,j,L1,l,S:
option remember:
k:=nops(Mod):
gu:=0:

if m>n then
 RETURN(0):
fi:

L1:=[]:
for l in IC do
L1:=[op(L1),[seq(l[i]-l[i+1],i=1..nops(l)-1)]]:
od:

S:=B:
for i from 1 to nops(IC) do
if evalb(m=IC[i][1]) then
if evalb(L1[i]=[]) then
RETURN(0):
else
S:=S union {L1[i]}:
fi:
fi:
od:


if m=n then
 RETURN(1):
fi:

if k=0 then 
for m1 from 1 to m do
if not member([],Yeled(m-m1,A union S)) then
gu:=gu+GP(m1,n-m,A,Mod,Yeled(m-m1, A union S),IC):

fi:
od:
RETURN(gu):
fi:


i:=m mod k:

for m1 from 1 to m do
if not member([],Yeled(m-m1,A union Mod[i+1] union S)) then

gu:=gu+GP(m1,n-m,A,Mod,Yeled(m-m1, S union A union Mod[i+1]),IC):

fi:
od:

gu:


end:

#Gxn(n,A,Mod,B,IC): the number of 
# partitions of n avoiding the partition patterns A globally, patterns in Mod locally (according to mod conditions),
#and patterns in B at the beginning, in addition to having no sub-partitions in IC (initial conditions)

Gxn:=proc(n,A,Mod,B,IC) local m:
option remember:
add(GP(m,n,A,Mod,B,IC),m=1..n):
end:

#GxnSeq(N,A,Mod,B,IC): the first N terms (including n=0) of the sequence of numbers of 
# partitions of n avoiding the partition patterns A globally, patterns in Mod locally (according to mod conditions),
#and patterns in B at the beginning, in addition to having no sub-partitions in IC (initial conditions)

GxnSeq:=proc(N,A,Mod,B,IC):
[1,seq(Gxn(n,A,Mod,B,IC),n=1..N)]:
end:


######################Search#########################

Findprod:=proc(n,A,Mod,B,IC) local LTS,q:
option remember:

LTS:=ListToSeries(GxnSeq(n,A,Mod,B,IC),q):

prodmake(LTS,q,n,list):

end:
#########################################################################################################################
#Search: input:
#N:we search for restricted partitions that that the product side matches the sum side up to at least N terms 
#(we then later verify they match to at least 100 terms)
#A: a set of sets (of the global conditions to avoid, we search through all the sets in A)
#Mod: the Mod condition to avoid 
#B: the set of beginning conditions to avoid (ususally we just take it to be {})
#IC: the list of initial conditions to avoid, for example, IC=[[1],[2]] indicates avoiding [1] and [2] as sub-partitions.
#S: the set of numbers allowed to appear in the list notation of the product side. We usually take this to be {1,0,-1,-2}. 

#output: S11 where
#S11: {[[A,ini,Mod,period,first 30 terms of the sequence],`Product side is:`, S, `mod `, P]} 
#(gathers all such partitions with product side period >= 4)


#local variables:
#Prod: {`Product side is:`, S, `mod `, P} (gathers all product sides with period >= 4)
#S12: {[[A,ini,Mod,period],`Product side is:`, S, `mod `, P]} (gathers all such partitions with product side period >= 10)
#Prod1: {`Product side is:`, S, `mod `, P} (gathers all product sides with period >= 10)
#S1: {[A,ini,Mod,period,first 30 terms of the sequence]}
#S2: {[first 30 terms of the sequence, period]}
###########################################################################################################################

Search:=proc(N,A,Mod,B,IC,S) local Prod,S11,Prod1,S12,A1,IC1,Modp,L,a,ini,i,k,s, m1,m2,m3,m4,co,Rec,S1,S2:
option remember:
S1:={}:
S2:={}:
A1:={}:


for a in A do 
if member([0], a) or member([0,0],a) or member([0,0,0],a) or member([0$4],a) or member([0$5],a) then
A1:=A1 union {a}:
fi:
od:
IC1:=CleanIC(IC):
IC1:=powerset(IC1):

k:=nops(Mod):
Modp:=[1$k]:
for i from 1 to k do 
Modp[i]:=powerset(Mod[i]):
od:

co:=0:
if k=0 then 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],[], PerFinder(L),GxnSeq(30,Rec[1],Rec[2],{},Rec[3])]}:
S2:=S2 union {[GxnSeq(30,Rec[1],[],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
fi:
od:
od:
fi:


co:=0:
if k=1 then
for m1 in Modp[1] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[m1],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L),GxnSeq(30,Rec[1],Rec[2],{},Rec[3])]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
fi:



co:=0:
if k=2 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[m1,m2],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2],PerFinder(L),GxnSeq(30,Rec[1],Rec[2],{},Rec[3])]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
fi:


co:=0:
if k=3 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[m1,m2,m3],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2],PerFinder(L),GxnSeq(30,Rec[1],Rec[2],{},Rec[3])]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
fi:

co:=0:
if k=4 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for m4 in Modp[4] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[m1,m2,m3,m4],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2],PerFinder(L),GxnSeq(30,Rec[1],Rec[2],{},Rec[3])]}: 
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
od:
fi:

#S11: {[[A,ini,Mod,period,first 30 terms of the sequence],`Product side is:`, S, `mod `, P]} 
#(gathers all such partitions with product side period >= 4)

#Prod: {`Product side is:`, S, `mod `, P} (gathers all product sides with period >= 4)


S11:={}:
Prod:={}:
for i from 1 to nops(S1) do
if S1[i][4]>=4 then
if not Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])=`Product side is not desired` then
S11:=S11 union {[S1[i],Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])]}:
Prod:=Prod union {Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])}:
fi:
fi:
od:

#S12: {[[A,ini,Mod,period,first 30 terms of the sequence],`Product side is:`, S, `mod `, P]} 
#(gathers all such partitions with product side period >= 10)

#Prod1: {`Product side is:`, S, `mod `, P} (gathers all product sides with period >= 10)


S12:={}:
Prod1:={}:
for i from 1 to nops(S1) do
if S1[i][4]>=10 then
if not Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])=`Product side is not desired` then
S12:=S12 union {[S1[i],Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])]}:
Prod1:=Prod1 union {Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])}:
fi:
fi:
od:


#Recall S11: {[[A,ini,Mod,first 30 terms of the sequence],`Product side is:`, S, `mod `, P]} 
#(gathers all such partitions with product side period >= 4)

for s in S11 do
print(s):
od:

end:







############################
#A variation of Search     #
#takes A and M as powersets#
############################

SearchAM:=proc(N,A,Mod,B,IC,S) local A1,IC1,Modp,L,a,ini,i,k, m1,m2,m3,m4,co,Rec,S1,S2,S11,Prod,S12,Prod1:
option remember:
S1:={}:
S2:={}:

A1:={}:

if A={} then
A1:={{}}:
fi:


for a in A do 
if member([0], a) or member([0,0],a) or member([0,0,0],a) or member([0$4],a) or member([0$5],a) then
A1:=A1 union {a}:
fi:
od:

IC1:=CleanIC(IC):


IC1:=powerset(IC1):

Modp:=Mod:
k:=nops(Mod):


co:=0:
if k=0 then 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[],ini):
L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then
S1:=S1 union {[Rec[1],Rec[3],[], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],[],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
fi:
od:
od:
fi:




co:=0:
if k=1 then
for m1 in Modp[1] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:

Rec:=Filter1(a,[m1],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
fi:



co:=0:
if k=2 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[m1,m2],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
fi:


co:=0:
if k=3 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter1(a,[m1,m2,m3],ini):  

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
fi:

co:=0:
if k=4 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for m4 in Modp[4] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi: 
Rec:=Filter1(a,[m1,m2,m3,m4],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:  #{[A, ini, Mod, period]}
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
od:
fi:

S11:={}:
Prod:={}:
print(`a and ini and Mod and period are:`):
for i from 1 to nops(S1) do
if S1[i][4]>=11 then 
if not Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])=`Product side is not desired` then
print(S1[i]):
S11:=S11 union {[S1[i],Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])]}:
Prod:=Prod union {Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])}:
fi:
fi:
od:



S12:={}:
Prod1:={}:
print(`a and ini and Mod and period are:`):
for i from 1 to nops(S1) do
if S1[i][4]>=11 then
if not Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])=`Product side is not desired` then
S12:=S12 union {[S1[i],Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])]}:
Prod1:=Prod1 union {Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])}:
fi:
fi:
od:

print(`the number of such partitions of n and period are:`):
for i from 1 to nops(S2) do
if S2[i][3]>=11 then
print(S2[i]):
fi:
od:

[S11,Prod]:


end:







#output: `Product side is:`, S, `mod `, P

Prodside:=proc(n,A,Mod,B,IC) local S,i,L,P,L1:

L:=Findprod(n,A,Mod,B,IC):

P:=PerFinder(L):

L1:=L[1..P]:

if not convert(L1,set) union {0,-1,-2,1}={0,-1,-2,1} then
RETURN(`Product side is not desired`):
fi:


if convert(L1,set) union {0,-1,-2,1}={-2,-1,0,1} then

if member(-2,L1) then
RETURN(`Product side is:`, L1):
fi:


if member(1,L1) then
RETURN(`Product side is:`, L1):
fi:

fi:

S:={}:

for i from 1 to P do
if L[i]=-1 then 
S:=S union {i}:
fi:
od:

RETURN(`Product side is:`, S, `mod `, P):

end:


#input a set A and a positive integer l, output the set of subsets of A whose length is no greater than l

powerL:=proc(A,l) local S1,S2,a:

S1:=powerset(A):

S2:={}:

for a in S1 do
if nops(a)<=l then
S2:=S2 union {a}:
fi:
od:

S2:

end:


powerLF:=proc(A,l) local S1,n,i,j,k:
option remember:

n:=nops(A):
S1:={{}}:

for i from 1 to n do
S1:=S1 union {{A[i]}}: 
od:

if l=2 then
for i from 1 to n do 
for j from i+1 to n do 
S1:=S1 union {{A[i],A[j]}}:
od:
od:
fi:

if l=3 then 

for i from 1 to n do 
for j from i+1 to n do 
S1:=S1 union {{A[i],A[j]}}:
od:
od:

for i from 1 to n do 
for j from i+1 to n do 
for k from j+1 to n do
S1:=S1 union {{A[i],A[j],A[k]}}:
od:
od:
od:
fi:


S1:

end:
#############################################
#all compositions with m parts, largest part=n (allowing 0 as a part) 

Comp1:=proc(m,n) local S,m1,p,L,p1:
option remember:

if n=0 then
RETURN({[0$m]}):
fi:

S:={}:

for m1 from 1 to m do 
for p in `union`(seq(Comp1(m-m1,n1),n1=0..n)) do
S:=S union {op(permute([op(p),n$m1]))}:
od:
od:


RETURN(S):
end:




#all compositions with <=m parts, largest part <=n (allowing 0 as a part)

Comp:=proc(m,n) local m1, n1,S:
option remember:
S:={}:
for m1 from 1 to m do
for n1 from 0 to n do
S:=S union Comp1(m1,n1):
od:
od:

S:
end:


####################################################

Filter1:=proc(A,Mod,IC) local a, A1, Mod1, Mod2, IC1,i,k,ModInter,ini,iniP:
option remember:


k:=nops(Mod):
Mod1:=Mod:

A1:=FilterS(A):
IC1:=FilterS(IC):
for i from 1 to k do
Mod1[i]:=FilterS(Mod[i]):
od:



#######if Mod has only one element, move it to A#######
if k=1 then
A1:=A1 union op(Mod1): 
Mod1:=[]:
k:=k-1:
fi:

if not (k=0 or k=1) then

###Getting rid of the elements of A in Mod########
for a in A1 do
for i from 1 to k do
if a in Mod1[i] then
Mod1[i]:=Mod1[i] minus {a}:
fi:
od:
od:


####Moving the common elements of Mod to A#######
ModInter:=`intersect`(seq(Mod1[i],i=1..k)):
A1:=A1 union ModInter:

for i from 1 to k do
Mod1[i]:=Mod1[i] minus ModInter:
od:


######If Mod include only {}, change Mod to []##########
Mod2:=`union`(seq(Mod1[i],i=1..k)):

if evalb(Mod2={}) then 
Mod1:=[]:
fi:

fi:


IC1:=convert(IC1,set):


#######Remove the members of IC that are redundent (because of A and Mod)######
for ini in IC1 do 
iniP:=PartoP(ini):
if member(iniP,A1) or member(iniP,Mod1) then 
IC1:=IC1 minus {ini}:
fi:
od:

IC1:=convert(IC1,list):
A1:=FilterS(A1):
RETURN([A1,Mod1,IC1]):

end:


#Input a partition, output the underlying partition pattern

PartoP:=proc(p) local L,i:

L:=[1$(nops(p)-1)]:
for i from 1 to nops(p)-1 do
L[i]:=p[i]-p[i+1]:
od:

L:
end:


########################
#A variation of Search #
#Uses fictitious zeros #
########################

SearchFic:=proc(N,A,Mod,B,numFic,S) local  A1,IC,IC1, Modp,L,a,ini,i,k, m1,m2,m3,co,Rec,S1,S2:
S1:={}:
S2:={}:

IC:=convert(Fic(numFic,A,Mod),list):

A1:=powerset(A):
IC1:=powerset(IC):

k:=nops(Mod):
Modp:=[1$k]:
for i from 1 to k do 
Modp[i]:=powerset(Mod[i]):
od:

co:=0:
if k=0 then 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi: 
Rec:=Filter(a,[],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],[], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],[],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
fi:
od:
od:
fi:


co:=0:
if k=1 then
for m1 in Modp[1] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter(a,[m1],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
fi:



co:=0:
if k=2 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter(a,[m1,m2],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
fi:


co:=0:
if k=3 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter(a,[m1,m2,m3],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
fi:

print(`a and ini and Mod and period are:`):
for i from 1 to nops(S1) do
if S1[i][4]>=12 then
print(S1[i]):
fi:
od:
print(`the number of such partitions of n and period are:`):
for i from 1 to nops(S2) do
if S2[i][3]>=12 then
print(S2[i]):
fi:
od:



end:






############
Search2:=proc(N,A,Mod,B,numFic,S) local  A1,IC,IC1, Modp,L,a,ini,i,k, m1,m2,m3,co,Rec,S1,S2:
S1:={}:
S2:={}:

IC:=convert(Fic(numFic,A,Mod),list):

A1:=powerset(A):
IC1:=powerset(IC):

k:=nops(Mod):
Modp:=[1$k]:
for i from 1 to k do 
Modp[i]:=powerset(Mod[i]):
od:

co:=0:
if k=0 then 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
if co>4000 then break: fi: 
Rec:=Filter(a,[],ini):
L:=Findprod(N,Rec[1],[],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then
S1:=S1 union {[Rec[1],Rec[3],[], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],[],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
od:
od:
fi:


co:=0:
if k=1 then
for m1 in Modp[1] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
if co>4000 then break: fi: 
Rec:=Filter(a,[m1],ini):
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then
S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
od:
od:
od:
fi:



co:=0:
if k=2 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
if co>4000 then break: fi: 
Rec:=Filter(a,[m1,m2],ini):
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then
S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
od:
od:
od:
od:
fi:


co:=0:
if k=3 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
if co>4000 then break: fi: 
Rec:=Filter(a,[m1,m2,m3],ini):
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then
S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
od:
od:
od:
od:
od:
fi:

print(`a and ini and Mod and period are:`):
for i from 1 to nops(S1) do
print(S1[i]):
od:

for i from 1 to nops(S2) do
print(S2[i]):
od:



end:






############
#m: number of fictitious 0s
#output: ICF

Fic:=proc(m,A,Mod) local a,M,ICF,p,k,k1,k2,i,Mod1,b,c:
ICF:={}:

if m=0 then
RETURN({}):
fi:


for a in A do

if not member(0,a) then 
k:=nops(PtoPar(a,0)):
ICF:=ICF union {[op(1..k-1,PtoPar(a,0))]}:
fi:

if numZero(a)<m then
p:=cutZero(a):
if not p=[] then
k1:=nops(PtoPar(p,0)):
ICF:=ICF union {[op(1..k1-1,PtoPar(p,0))]}:
fi:
fi:

od:


k2:=nops(Mod):
Mod1:=[{}$k2]:

for i from 1 to k2 do
for b in Mod[i] do

if not member(0,b) and (add(b) mod k2 =i-1) then 
k:=nops(PtoPar(b,0)):
Mod1[i]:=Mod1[i] union {[op(1..k-1,PtoPar(b,0))]}:
fi:

if numZero(b)<m and (add(b) mod k2 =i-1) then
p:=cutZero(b):
if not p=[] then
k1:=nops(PtoPar(p,0)):
Mod1[i]:=Mod1[i] union {[op(1..k1-1,PtoPar(p,0))]}:
fi:
fi:

od:
od:

`union`(seq(c,c in Mod1),ICF):


end:


#cut the zeros at the end of the list p
cutZero:=proc(p) local k:
k:=numZero(p):

[op(1..nops(p)-k,p)]:
end:




#CleanIC: clean the initial condition, get rid of elements that contain 0, 
#Also get rid of elements that are not actually a partition (order not correct)
#Input IC looks like [[1],[1,1],[2,1],[0],[2,3]], the output should be [[1], [1, 1], [2, 1]]

CleanIC:=proc(IC) local S,a:

S:=convert(IC,set):
for a in S do
if member(0,a) then
S:=S minus {a}:
fi:

if not IsPar(a) then 
S:=S minus {a}:
fi:
od:

convert(S,list):
end:

##
IsPar:=proc(p) local i:

for i from 1 to nops(p)-1 do
if p[i]<p[i+1] then
RETURN(false):
fi:
od:

true:
end:


#FilterS(A): input a set A of lists, filter out the sublists 
#For example, if A={[1,1,2],[1,2]} then the output will be {[1,1,2]}.



FilterS:=proc(A) local a,b,S:
S:=convert(A,set):

for a in A do
for b in A do
if not a=b then
if IsSubL(a,b) then
S:=S minus {a}:
fi:
fi:

od:
od:

if type(A,list) then
RETURN(convert(S,list)):
fi:

S:

end:


#IsSubL(L1,L2): takes two lists L1 and L2 and outputs true if L2 is a sublist of L1, false if not.

IsSubL:=proc(L1,L2) local i,j, n, k:

n:=nops(L1): 
k:=nops(L2):

for i from 1 to n-k+1 do 

for j from 1 to k do

if not L1[i+j-1]=L2[j] then
break:
fi:

if j=k then 
RETURN(true):
fi:

od:

od:

RETURN(false):
end:


#Poch(a,q,n) : the Pochhammer symbol expression

Poch:=proc(a,q,n):

product(1-a*q^i,i=0..n-1):

end:



ListToSeries:=proc(L,q)

return add(L[i]*q^(i-1),i=1..nops(L)):

end:


# IsPer(L,k):
# Inputs a list L
# and a possible period k
# Returns true if L is periodic
# with period k
#
# Try: IsPer([1,2,3,4,1,2,3],4);

IsPer:=proc(L,k) local i:
for i from k+1 to nops(L) do
	if L[i]<>L[i-k] then
		return false:
	fi:
od:
return true:
end:

# PerFinder(L):
# Finds the smallest period of a list L
# (note that if L is not periodic, it will return nops(L) )
#
# Try: PerFinder([1,2,1,2,1,2,1,2]);

PerFinder:=proc(L) local i:
for i from 1 to nops(L) do
	if IsPer(L,i) then
		return i:
	fi:
od:
return nops(L):
end:



#######################
#Nandi type conditions#
#######################

#takes A and M as a powerset, with Nandi-type conditions

SearchN:=proc(N,A,Mod,B,IC,S) local A1,IC1,Modp,L,a,ini,i,k, m1,m2,m3,m4,co,Rec,S1,S2,S11,Prod,S12,Prod1:
option remember:
S1:={}:
S2:={}:
A1:={}:

if A={} then
A1:={{}}:
fi:


for a in A do 
#if member([0], a) or member([0,0],a) or member([0,0,0],a) or member([0$4],a) or member([0$5],a) then
A1:=A1 union {a}:
#fi:
od:

IC1:=IC:
#IC1:=CleanIC(IC):


#IC1:=powerset(IC1):

Modp:=Mod:
k:=nops(Mod):


co:=0:
if k=0 then 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter2(a,[],ini):
L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then
S1:=S1 union {[Rec[1],Rec[3],[], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],[],{},Rec[3]), `period`,PerFinder(L)]}:
fi:
fi:
od:
od:
fi:




co:=0:
if k=1 then
for m1 in Modp[1] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:

Rec:=Filter2(a,[m1],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
print(L):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
fi:



co:=0:
if k=2 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter2(a,[m1,m2],ini):
L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
fi:


co:=0:
if k=3 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi:
Rec:=Filter2(a,[m1,m2,m3],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
fi:

co:=0:
if k=4 then
for m1 in Modp[1] do
for m2 in Modp[2] do 
for m3 in Modp[3] do 
for m4 in Modp[4] do
for a in A1 do
for ini in IC1 do
co:=co+1:
if co mod 1000 =0 then print(co): fi: 
Rec:=Filter2(a,[m1,m2,m3,m4],ini):

L:=Findprod(30,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S then
L:=Findprod(N,Rec[1],Rec[2],B,Rec[3]):
if {op(L),op(S)}=S and PerFinder(L)<=nops(L)/2 then

S1:=S1 union {[Rec[1],Rec[3],Rec[2], PerFinder(L)]}:
S2:=S2 union {[GxnSeq(30,Rec[1],Rec[2],{},Rec[3]), `period`,PerFinder(L)]}:

fi:
fi:
od:
od:
od:
od:
od:
od:
fi:



S11:={}:
Prod:={}:
print(`a and ini and Mod and period are:`):
for i from 1 to nops(S1) do
if S1[i][4]>=11 then 
if not Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])=`Product side is not desired` then
print(S1[i]):
S11:=S11 union {[S1[i],Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])]}:
Prod:=Prod union {Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])}:
fi:
fi:
od:



S12:={}:
Prod1:={}:
print(`a and ini and Mod and period are:`):
for i from 1 to nops(S1) do
if S1[i][4]>=11 then
if not Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])=`Product side is not desired` then
S12:=S12 union {[S1[i],Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])]}:
Prod1:=Prod1 union {Prodside(100,S1[i][1],S1[i][3],{},S1[i][2])}:
fi:
fi:
od:

print(`the number of such partitions of n and period are:`):
for i from 1 to nops(S2) do
if S2[i][3]>=11 then
print(S2[i]):
fi:
od:

[S11,Prod,S12,Prod1]:


end:

####################################################
Filter2:=proc(A,Mod,IC) local a, A1, Mod1, Mod2, IC1,i,k,ModInter,ini,iniP:
option remember:


k:=nops(Mod):
Mod1:=Mod:

A1:=FilterS(A):
IC1:=FilterS(IC):
for i from 1 to k do
Mod1[i]:=Mod[i]:
od:

#######if Mod has only one element, move it to A#######
if k=1 then
A1:=A1 union op(Mod1): 
Mod1:=[]:
k:=k-1:
fi:

if not (k=0 or k=1) then

###Getting rid of the elements of A in Mod########
for a in A1 do
for i from 1 to k do
if a in Mod1[i] then
Mod1[i]:=Mod1[i] minus {a}:
fi:
od:
od:


####Moving the common elements of Mod to A#######
ModInter:=`intersect`(seq(Mod1[i],i=1..k)):
A1:=A1 union ModInter:

for i from 1 to k do
Mod1[i]:=Mod1[i] minus ModInter:
od:


######If Mod include only {}, change Mod to []##########
Mod2:=`union`(seq(Mod1[i],i=1..k)):

if evalb(Mod2={}) then 
Mod1:=[]:
fi:

fi:


IC1:=convert(IC1,set):


#######Remove the members of IC that are redundent (because of A and Mod)######
for ini in IC1 do 
iniP:=PartoP(ini):
if member(iniP,A1) or member(iniP,Mod1) then 
IC1:=IC1 minus {ini}:
fi:
od:

IC1:=convert(IC1,list):


A1:=FilterS(A1):


RETURN([A1,Mod1,IC1]):



end:





################

#Star(L,N): input a list with Star (stand for *)  and an integer N, output a set of lists incroperating the * condition in Nandi-type partitions
#For example, Star([3,S,2,3,0], 30,1) outputs {[3,3,0],[3,2,3,0],[3,2,2,3,0],[3,2,2,2,3,0]}, 
#Star([3,Star,2,3,0], 30,2) outputs {[3,3,0],[3,2,2,3,0]}

Star:=proc(L,N,t) local k,i,S1,m,l,j:

k:=nops(L):
S1:={}:

for i from 1 to k do 
if L[i]=Star then 
for m from 0 to N do

l:=[op(1..i-1,L),L[i+1]$m*t,op(i+2..k,L)]:

if add(add([op(j..nops(l),l)]),j=1..nops(l))<N then 
S1:=S1 union {l}:
else 
break:
fi:


od:
fi:
od:
S1:
end:




#Comb(P,S): input a powerset P and a set S, create a new powerset, with S unioned with each set of P. 

Comb:=proc(P,S) local s,P1:

P1:={}:

for s in P do

P1:=P1 union {s union S}

od:

P1:
end: