#Version of March 18, 2019

print(`This is rPar.txt, to study and enumerate r-partitions written by Mingjia Yang`):
print(`It accompanies the article: `):
print(`Counting  r-Partitions by  Mingjia Yang`):
print(`For a list of the MAIN procedures type: `):
print(`Help();`):
print(`For a specific procedure, type Help(ProcedureName); For example, Help(Pr); `):
print(`For a list of the supporting procedures, type: Help1();`):
print(`For a specific procedure, type Help(ProcerureName); For example, Help(Ar); `):
with(ListTools):


Help1:=proc():

if args=NULL then
 print(`The supporting procedures are:  Ar, Ch, Ftrq, FtrqN,GenPr,  GenRP1, GenRPr, NAr, RL, RP,   RPr, RPr1`):
 print(` `):
else
Help(args):
fi:
end:


Help:=proc():

if args=NULL then
 print(`The main procedures are: checkF, F, Fr, GuessPol, NPr, NPrSeq, Pr, qGuessPol, Poch,SubL,SubLk,Ex,G,GP` ):
 print(` `):


elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(L,n,s0): guesses a polynomial in n for`):
print(`the list L, such that L[i]=P[i] for i>=s0 for example, try: `):
print(`GuessPol([seq(i,i=1..10)],n,1);`):

elif nops([args])=1 and op(1,[args])=GP then
print(`GP(M,N) input integers M and N`):
print(`output the Gaussian polynomial (M,N)_q `):
print(`Try: `):
print(`GP(2,1)`):

elif nops([args])=1 and op(1,[args])=G then
print(`G(M,N) input integers M and N`):
print(`output GP(2N+M-1,N) `):
print(`Try: `):
print(`G(2,3)`):

elif nops([args])=1 and op(1,[args])=Pr then
print(`Pr(n,r): takes two integers n and r and outputs the set of r-partitions (a[i+1]<=a[i]+r) of n `):
print(`Try: `):
print(`Pr(5,2);`):


elif nops([args])=1 and op(1,[args])=NPrSeq then
print(`NPrSeq(N,r): the the first N+1 terms of the sequence NPr(n,r). Try`):
print(`NPrSeq(20,1);`):

elif nops([args])=1 and op(1,[args])=Ar then
print(`Ar(n,m,r): takes three integers n, m and r, outputs the set of r-partitions of n with the first part no greater than m. `):
print(`Try: `):
print(`Ar(10,4,1);`):


elif nops([args])=1 and op(1,[args])=NPr then
print(`NPr(n,r): takes two integers n and r and outputs the number of r-partitions of n (that is, NPr is the numeric version of Pr)`):
print(`Try: `):
print(`NPr(5,2);`):


elif nops([args])=1 and op(1,[args])=NAr then
print(`NAr(n,m,r): the numeric version of Ar(n,m,r) `):
print(`Try: `):
print(`NAr(10,4,1);`):

elif nops([args])=1 and op(1,[args])=RP then
print(`RP(n,M,N): the set of restrained partitions with the first part no greater than M and the number of parts no greater than N.`):
print(`Try: `):
print(`RP(8,5,4);`):

elif nops([args])=1 and op(1,[args])=RL then
print(`RL(S,N): takes a set of lists and filters it so that only the lists with length no greater than N remains. `):
print(`Try: `):
print(`RL({[1,2,3,4],[1,2],[4]},2);`):

elif nops([args])=1 and op(1,[args])=GenRP then
print(`GenRP(k,M,N,q): (using RP) outputs the first k terms of the generating function for restrained partitions with the first part <=M and the number of parts<=N.`):
print(`Try: `):
print(`GenRP(20,5,4,q);`):


elif nops([args])=1 and op(1,[args])=GenRP1 then
print(`GenRP1(k,M,N,q): using the known formula product(1-q^(M+i),i=1..N)/product(1-q^j,j=1..N), compute the first k terms of the generating function for restrained partitions with the first part <=M and the number of parts<=N.`):
print(`Try: `):
print(`GenRP1(20,5,4,q);`):

elif nops([args])=1 and op(1,[args])=RPr then
print(`RPr(n,r,M,N): outputs the set of r-partitions of n with the first part no greater than M and the number of parts no greater than N.`):
print(`Try: `):
print(`RPr(10,1,5,4)`):


elif nops([args])=1 and op(1,[args])=GenRPr then
print(`GenRPr(k,r,M,N,q): (using RRP) outputs the first k terms of the generating function for restrained r-partitions with the first part <=M and the number of parts<=N.`):
print(`Try: `):
print(`GenRPr(20,1,5,4,q)`):


elif nops([args])=1 and op(1,[args])=F then
print(`F(M,N,r,q): The generating function for all r-partitions with EXACTLY N parts and first part EXACTLY M. Try:`):
print(`F(4,3,1,q);`):

elif nops([args])=1 and op(1,[args])=Fr then
print(`Fr(M,N,r,q): The generating function for all r-partitions with<= M parts and the first part<= N. Try:`):
print(`Fr(4,3,1,q);`):

elif nops([args])=1 and op(1,[args])=Ch then
print(`Ch(M,N,r):  Checks Fr(M,N,r,q) against the brute-force version GenPr. Try:`):
print(`Ch(4,3,1);`):

elif nops([args])=1 and op(1,[args])=RPr1 then
print(`RPr1(M,N,r): the set of r-partitions with at most N parts and the first part <=M` ):
print(`RPr1(2,3,0);`):

elif nops([args])=1 and op(1,[args])=GenPr then
print(`GenPr(N,r,q): the the first N+1 terms of the generating function of NPr(n,r)`):
print(`GenPr(20,1,q);`):

elif nops([args])=1 and op(1,[args])=GenPrSeq then
print(`GenPrSeq(r,N): the first N terms of the sequence of -(r-1)-partitions. `):


elif nops([args])=1 and op(1,[args])=checkF then
print(`checkF(M,N,r): check if the conjectured formula corresponding to F(M,N,r,1) is true  `):
print(`Try: checkF(M,N,r)`):


elif nops([args])=1 and op(1,[args])=qGuessPol then
print(`qGuessPol(L,n,q,s0): guesses a polynomial in q^n for the list L, such that L[i]=P[i] for i>=s0`):
print(`Try: qGuessPol([seq(F(M1,2,-1,q),M1=1..20)],M,q,1)`):


elif nops([args])=1 and op(1,[args])=Ftrq then
print(`Ftrq(M,N,q): outputs the formula for Ftr(M,N,-1,q) guessed using qGuessPol, N<=5`):
print(`Try: Ftrq(M,4,q)`):

elif nops([args])=1 and op(1,[args])=FtrqN then
print(`Ftrq(M,N,q): outputs the numerator of Ftrq (N<=5)`):
print(`Try: FtrqN(M,4,q)`):

elif nops([args])=1 and op(1,[args])=Poch then
print(`Poch(a,q,n): the Pochhammer symbol `):
print(`Try: `):
print(`Poch(a,q,3)`):

elif nops([args])=1 and op(1,[args])=SubL then
print(`SubL(L): inputs a list and outputs a list that is the consecutive difference of elements in L. `):
print(`Try:  `):
print(`SubL([1, 2, 4, 7, 11, 16, 22, 29, 37])`):


elif nops([args])=1 and op(1,[args])=SubLk then
print(`SubLk(L,k): SubL applied k times to L.  `):
print(`Try:`):
print(`SubL([1, 2, 4, 7, 11, 16, 22, 29, 37],2)`):


elif nops([args])=1 and op(1,[args])=Ex then
print(`Ex(L,i):  input a list of polynomials, extract the coeffient of the ith term from the last (including the last) of each of these polynomials and return the list of these coefficients `):
print(`Try:`):
print(`Ex([seq(F(2,j,-1,q),j=5..20)],6);`):



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




#GuessPol(L,n,s0): guesses a polynomial in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol([(seq(i,i=1..10)],n,1);
GuessPol:=proc(L,n,s0) local d,gu:

for d from 0 to nops(L)-s0-3 do
 gu:=GuessPol1(L,d,n,s0):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:

#GuessPol1(L,d,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol1([(seq(i,i=1..10),1,n,1);
GuessPol1:=proc(L,d,n,s0) local P,i,a,eq,var:
if d>=nops(L)-s0-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=s0..s0+d+3)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:



#Pr(n,r) takes two integers n and r and outputs the set of r-partitions of n

Pr:=proc(n,r) local m:

`union`(seq(Ar(n,m,r),m=1..n)):

end:




#Ar(n,m,r) takes three integers m and n and outputs the set of r-partitions of n with 
#the first part no greater than m. 

Ar:=proc(n,m,r) local m1,S,i:
option remember:

S:={}:

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

if m>n then
RETURN({}):
fi:

for m1 from 1 to m-r do
S:=S union {seq([m,op(Ar(n-m,m1,r)[i])],i=1..nops(Ar(n-m,m1,r)))} :
od:

S:
end:


#NPr(n,r) takes two integers n and r and outputs the number of r-partitions of n

NPr:=proc(n,r) local m:
add(NAr(n,m,r),m=1..n):
end:


#NPrSeq(N,r): the first N+1 terms of the generating function of NPr(n,r)
NPrSeq:=proc(N,r) local n1:

[1,seq(NPr(n1,r),n1=1..N)]:

end: 

#GenPr(N,r,q): the first N+1 terms of the generating function of NPr(n,r)
GenPr:=proc(N,r,q):
1+add(NPr(n,r)*q^n,n=0..N):
end: 


#NAr(n,m,r) takes three integers m and n and outputs the set of r-partitions of n 
#with the the first part being m. 

NAr:=proc(n,m,r) local m1,S,i:
option remember:

S:=0:

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

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

for m1 from 1 to m-r do
S:=S + NAr(n-m,m1,r) :
od:

S:
end:




#RP(n,M,N): the set of resticted partitions with the first part no greater than M 
#and the number of parts no greater than N.


RP:=proc(n,M,N) local m,S,i: 
S:={}:
for m from 1 to M do
S:=S union Ar(n,m,0):
od:

RL(S,N):
end:


#RL(S,N) takes a set of lists and filters it so that only the lists with length no greater than N remains. 

RL:=proc(S,N) local s,S1:
S1:={}:
for s in S do
if nops(s)<=N then
S1:=S1 union {s}:
fi:
od:

S1:
end:



#GenRP: the first k terms of the generating function for restrained partitions 
#with the the first part <=M and the number of parts<=N. 

GenRP:=proc(k,M,N,q) local n:

1+add(nops(RP(n,M,N))*q^n,n=0..k):

end:

#GenRP1: use the known formula to compute the same thing as GenP

GenRP1:=proc(k,M,N,q) local i,j:

taylor((product(1-q^(M+i),i=1..N)/product(1-q^j,j=1..N)),q=0,k+1):

end:



#RPr(n,r,M,N): outputs the set of r-partitions of n with the first part no greater than M 
#and the number of parts no greater than N.

RPr:=proc(n,r,M,N) local m,S,i: 
S:={}:
for m from 1 to M do
S:=S union Ar(n,m,r):
od:

RL(S,N):
end:

#GenRPr: the first k terms of the generating function for restrained r-partitions 
#with the the first part <=M and the number of parts<=N. 

GenRPr:=proc(k,r,M,N,q) local n:

1+add(nops(RPr(n,r,M,N))*q^n,n=0..k):

end:




#F(M,N,r,q): The generating function for all r-partitions with EXACTLY N parts and FIRST part EXACTLY M. Try:
#F(4,3,1,q);
F:=proc(M,N,r,q) local M1:
option remember:

if N=1 then
q^M:
else
 expand(q^M*add(F(M1,N-1,r,q),M1=1..M-r)):
fi:
end:



#FtrNew(M,N,r,q): The generating function for all r-partitions with EXACTLY N parts and FIRST part EXACTLY N. Try:
#FtrNew(4,3,1,q);
FtrNew:=proc(M,N,r,q) 
option remember:
if r>0 then
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:

if N=1 then
RETURN(q^M):
else
 if M<1 then
  RETURN(0):
else
RETURN( expand( q*FtrNew(M-r,N,r,q)+ q^M*FtrNew(M-r,N-1,r,q) )):
fi:
fi:
end:

#Fr(M,N,r,q): The generating function for all r-partitions with at most M parts and the first part at most N. Try:
#Fr(4,3,1,q);
Fr:=proc(M,N,r,q) local M1,N1:
1+add(add(F(M1,N1,r,q),M1=1..M),N1=1..N):
end:

#Ch(M,N,r):Checks Fr(M,N,r,q) against the brute-force version
Ch:=proc(M,N,r) local q,gu; gu := Fr(M, N, r, q); 
evalb(GenRPr(degree(gu, q), r, M, N, q) - gu=0): 
end:



#the set of r-partitions with at most N parts and the first part part <=M 
RPr1:=proc(M,N,r) local n,i,S:
S:={}:

n:=M*N+r*N*(N-1)/2:

for i from 0 to n do
S:=S union RPr(i,r,M,N):
od:

S:
end:




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


#check if the conjectured formula corresponding to F(M,N,r,1) is true 


checkF:=proc(M,N,r) local i,j,M1,c: 

if N=1 then
RETURN(true):
fi:

c:=simplify((M-r)*product(i,i=M-r*N+1..M+(1-r)*N-2)/product(j,j=1..N-1)-sum((M1-r)*product(i,i=M1-r*(N-1)+1..M1+(1-r)*(N-1)-2)/product(j,j=1..N-2),M1=1..M-r)):

if c=0 then
RETURN(true):
fi:

false:

end:



#qGuessPol(L,n,q,s0): guesses a polynomial in q^n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#qGuessPol([(seq(i,i=1..10)],n,q,1);
qGuessPol:=proc(L,n,q,s0) local d,gu:

for d from 0 to nops(L)-s0-3 do
 gu:=qGuessPol1(L,d,n,s0):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:

#qGuessPol1(L,d,n,s0): guesses a polynomial of degree d in q^n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#qGuessPol1([seq(i,i=1..10)],1,n,q,1);
qGuessPol1:=proc(L,d,n,s0) local P,i,a,eq,var:
if d>=nops(L)-s0-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*(q^n)^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=s0..s0+d+3)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:


RecVerify:=proc(M,N) local c:

c:=simplify(qGuessPol([seq(F(M1,N,-1,q),M1=1..20)],M,q,1)-q*qGuessPol([seq(F(M1,N,-1,q),M1=1..20)],M-1,q,1)-q^M*qGuessPol([seq(F(M1,N-1,-1,q),M1=1..20)],M+1,q,1)):

if c=0 then 
RETURN(true):
fi:
c:
end:



#outputs the formula guessed using qGuessPol
Ftrq:=proc(M,N,q):

if N=1 then 
RETURN(q^M):
fi:

if N=2 then
RETURN(q^(M+1)*(q^(M+1)-1)/(q-1)):
fi:

if N=3 then
RETURN(q^(M+2)*(q^(M+3)+q^2-q-1)*(q^(M+1)-1)/(q-1)^2/(q+1)):
fi:

if N=4 then 
RETURN(q^(M+3)*(q^(2*M+8)+q^(M+7)-q^(M+5)-q^(M+4)-q^(M+3)+q^6-2*q^4-q^3+2*q+1)*(q^(M+1)-1)/(q-1)^3/(q+1)/(q^2+q+1)):
fi:

if N=5 then
RETURN(q^(M+4)*(q^(M+5)+q^4-q-1)*(q^(2*M+10)-q^(M+4)-q^(M+3)+q^8-q^5-2*q^4+2*q+1)*(q^(M+1)-1)/(q-1)^4/(q+1)^2/(q^2+1)/(q^2+q+1)):
fi:

end:


#the numerators for F(M,N,r,q)
FN:=proc(M,N,r,q):

expand(simplify(numer(F(M,N,r,q)))):
end:


GP:=proc(a,b) local i, j:

product(1-q^i,i=a-b+1..a)/product(1-q^j,j=1..b):
end:


G:=proc(M,N):

expand(simplify(GP(2*N+M-1,N))):

end: 

FtrqN:=proc(M,N,q):

expand(numer(Ftrq(M,N,q))):

end:





GuessB:=proc(M,N,r) local A,B:

A:=(M-r)*(M+(1-r)*N-2)!/(N-1)!/(M-r*N)!:

B:=binomial(M+(1-r)*N/(-r)-2,N*(-r)-1)-binomial(M+(1-r)*N/(-r)-2,N/(-r)+r-1):

[A,B]:
end:




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

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

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

end:


#SubL(L) inputs a list and outputs a list that is the consecutive difference of elements in L. 
SubL:=proc(L) local L1,i:

L1:=[0$(nops(L)-1)]:

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

L1:

end:

#SubLk(L,k) SubL applied k times to L. 

SubLk:=proc(L,k) local L1,i,L2: 
option remember:

if k=0 then
RETURN(L):
fi:

L1:=[0$(nops(L)-1)]:

L1:=SubL(L):
L2:=SubLk(L1,k-1):

RETURN(L2): 
end:


#Ex(L,i): input a list of polynomials, extract the coeffient of the ith term from the last of each of these polynomials
#and return the list of the ith terms. (the polynomials in L have to contain more than one term)

Ex:=proc(L,i) local L1,p:

L1:=[]:

for p in L do 
L1:=[op(L1),subs(q=1,convert(p,list)[-i])]:
od:
L1:

end: