######################################################################
##DMP: Save this file as  DMP                                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read DMP<Enter>                                    #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Jan. 18, 2012
#This version: Jan. 25, 2012 (adding RandDMP)
print(`In LOVING MEMORY of my MENTOR and GURU `):
print(` Herbert Saul WILF (June 13, 1931-Jan. 7, 2012)`):
print(``): 
print(`Created:  Jan. 12, 2012`):
print(`This version:  Jan. 25, 2012 [adding procedure RandDMP(n)] `):
print(` This is DMP `):
print(`A Maple package to study Distinct Multiplicity Partitions`):
print(`inspired by Problem #6 in Herb Wilf's on-line article `):
print(` "Some unsolved problems" `):
print(` Herb Wilf's article is available from: `):
print(`http://www.math.upenn.edu/%7Ewilf/website/UnsolvedProblems.pdf `):
print(``):
print(`The present package accompanies the article `):
print(`Using GENERATINGFUNCTIONOLOGY to Enumerate`):
print(`Distinct-Multiplicity Partitions`):
print(`by  Doron Zeilberger`):
print(`that is available from:`):
print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/dmp.html`):
print(`The url of the present package is:`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/DMP`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
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(`---------------------------------------------------`):
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 story procedures type ezraStories();,`):
print(`  for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------------------`):
with(combinat):

ezraStories:=proc()

if args=NULL then
 print(` The story  procedures are: `):
 print(` SipurChange, SipurDMP, SipurSeq, SipurSeqPC, `):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` Aitken, qnmS, RandDMP1, RO, LoadedDie `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` DMPS, GFmxq, GFmq, GFSq, GFSxq, pn, pnSeq, qn,qnPC,`):
print(` qnSeqPC, , RandDMP, SeqSN,  `):


elif nargs=1 and args[1]=Aitken then
print(`Aitken(L): The conjectured limit of the sequence L`):
print(`from its three last elements, using Aitken's fromula`):
print(`For example, try: Aitken([3.1,3.01,3.001]);`):

elif nops([args])=1 and op(1,[args])=DMPS then
print(`DMPS(S,n): the set of distinct-multiplicity partitions of`):
print(`n using the members of S using generating functions`):
print(`try DMPS([1,3,4],10);`):


elif nops([args])=1 and op(1,[args])=GFmq then
print(`GFmq(m,q): The full generating function for integer partitions`):
print(`using {1, ...,m} in terms of q`):
print(`The coeff. of q^n in the Maclaurin expansion of the output`):
print(`gives you the number of partitions of n into distinct multiplicities`):
print(`whose largest part is <=m`):
print(`For example, try:`):
print(`GFmq(4,q);`):



elif nops([args])=1 and op(1,[args])=GFmxq then
print(`GFmxq(m,x,q): The full generating function for integer partitions`):
print(`into distinct multiplicities`):
print(`using {1, ...,m} in terms of x[1], ..., x[m] and q`):
print(`For example, try:`):
print(`GFmxq(4,x,q);`):

elif nops([args])=1 and op(1,[args])=GFSq then
print(`GFSq(S,q): The full generating function for integer partitions`):
print(`using the parts in the list (or set) S, in terms of q`):
print(`The coeff. of q^n in the Maclaurin expansion of the output`):
print(`gives you the number of partitions of n into distinct multiplicities`):
print(`whose largest part is <=m`):
print(`For example, try:`):
print(`GFSq([1,5,10,25],q);`):

elif nops([args])=1 and op(1,[args])=GFSxq then
print(`GFSxq(S,x,q): The full generating function for integer partitions`):
print(`using the parts in the list (or set) S, in terms of q`):
print(`The coeff. of q^n in the Maclaurin expansion of the output`):
print(`gives you the weight-enumerator, in terms of the indexed variables`):
print(` x[i] `):
print(` of partitions of n into distinct multiplicities`):
print(`where the parts are drawn from the set S`):
print(`For example, try:`):
print(`GFSxq([1,5,10,25],x,q);`):

elif nops([args])=1 and op(1,[args])=pn then
print(`pn(n): the number of partitions of n`):
print(`For example, try:`):
print(` pn(20);`):

elif nops([args])=1 and op(1,[args])=pnSeq then
print(`pnSeq(n): the first N members of the sequence`):
print(`"number of multiplicity-free partitions of n"`):
print(`In other words OEIS sequence A098859. `):
print(`For example, try:`):
print(` pnSeq(30);`):

elif nops([args])=1 and op(1,[args])=qn then
print(`qn(n): the number of multiplicity-free partitions of n`):
print(`For example, try:`):
print(` qn(20);`):

elif nops([args])=1 and op(1,[args])=qnPC then
print(`qnPC(n): the pre-computed number of multiplicity-free partitions of n`):
print(`for n>220 we used Mciej Ireneusz Wilczynsk's table given at OEIS`):
print(`A098859 `):
print(`For example, try:`):
print(` qnPC(20);`):


elif nops([args])=1 and op(1,[args])=qnmS then
print(`qnmS(n,m,S): the number of partitions of n with largest part <=m`):
print(`with distinct multiplicities, none of which belongs to the set S`):
print(`For example, try:`):
print(`qnmS(10,10,{});`):

elif nops([args])=1 and op(1,[args])=qnSeq then
print(`qnSeq(n): the first N members of the sequence`):
print(`"number of multiplicity-free partitions of n"`):
print(`In other words OEIS sequence A098859. `):
print(`For example, try:`):
print(` qnSeq(30);`):

elif nops([args])=1 and op(1,[args])=qnSeqPC then
print(`qnSeqPC(n): the first N members of the sequence`):
print(`"number of multiplicity-free partitions of n"`):
print(`In other words OEIS sequence A098859. `):
print(`using pre-computed values`):
print(`For example, try:`):
print(` qnSeqPC(30);`):


elif nops([args])=1 and op(1,[args])=RandDMP then
print(`RandDMP(n): A uniformly-at-random distinct-multiplicity`):
print(`partition of n. It uses the methodology in `):
print(`the Nijenhuis-Wilf classic book Combinatorial Algorithms`):
print(`For example, try:`):
print(`RandDMP(10);`):

elif nops([args])=1 and op(1,[args])=RandDMP1 then
print(`RandDMP1(n,m,S): A uniformly-at-random distinct-multiplicity`):
print(`partition of n with largest part <=m`):
print(`none of which belongs to the set S`):
print(`For example, try:`):
print(`RandDMP1(10,10,{});`):

elif nops([args])=1 and op(1,[args])=RO then
print(`RO(F,m,r,t,q): applies the raising operator to`):
print(`F(q,t[1], ..., t[r]) that goes from`):
print(`the weight-enumerator of distinct-multiplicity partitions`):
print(`into exactly r parts using {1, ..., m-1} to those`):
print(`with exactly r+1 parts using {1, ..., m}, t is the`):
print(`symbol that carries the indexed variable for the parts with`):
print(`distinct multiplicities, q carries the number being partitioned`):
print(`For example, try`):
print(` RO(q^7*t[1]^7,2,1,t,q); `):

elif nops([args])=1 and op(1,[args])=SeqmN then
print(`SeqmN(m,N): the first N terms of the sequence`):
print(`"partitions of n into parts  <=m with distinct multiplicities"`):
print(`For example, try:`):
print(`SeqmN(4,40);`):

elif nops([args])=1 and op(1,[args])=SeqSN then
print(`SeqSN(S,N): the first N terms of the sequence`):
print(` "partitions of n into parts  drawn from the list S `):
print(` with distinct multiplicities" `):
print(`For example, try:`):
print(`SeqmN([1,5,10,25],40);`):

elif nops([args])=1 and op(1,[args])=SipurChange then 
print(`SipurChange(S,q,L): inputs a list (or set) of positive integers`):
print(` and outputs a webbook listing the (ordinary) generating function`):
print(`for the sequence "number of partitions of n`):
print(`with distinct multiplicity using the elements of S as parts`):
print(`It also lists the first L terms`):
print(`of each sequence`):
print(`For the sake of comparison it also does the (trivial) analogous`):
print(`thing when there is no restiction on the multipilicity`):
print(`For example, try:`):
print(`SipurChange([1,5,10,25],q,50);`):

elif nops([args])=1 and op(1,[args])=SipurDMP then 
print(`SipurDMP(M,q,L): inputs an integer N and outputs`):
print(`a webbook listing the (ordinary) generating functions`):
print(`for the sequence "number of partitions of n`):
print(`with distinct multiplicity using {1, ...,m}", for`):
print(`m from 1 to M. It also lists the first L terms`):
print(`of each sequence`):
print(`For example, try:`):
print(`SipurDMP(5,q,50);`):

elif nops([args])=1 and op(1,[args])=SipurSeq then 
print(`SipurSeq(L): inputs an integer L `):
print(` and outputs the first L terms of the`):
print(`sequence "number of partitions of n with distinct multiplicity"`):
print(`It also lists the first L terms`):
print(`of the classical sequence "number of partitions"`):
print(`and tries to find the asymptotic trend`):
print(`SipurSeq(50);`):

elif nops([args])=1 and op(1,[args])=SipurSeqPC then 
print(`SipurSeqPC(L): Like SipurSeq(L) but using precomputed values`):
print(`by Maciej Ireneusz Wilczynsk .`):
print(`L must be <=508`):
print(`Try: `):
print(`SipurSeqPC(508);`):

elif nops([args])=1 and op(1,[args])=LoadedDie then 
print(`LoadedDie(n,m,S): picks the multiplicity of m in a random`):
print(`distinct-multiplicity partition of n with largest part`):
print(`<=m and avoiding multiplicities in S. For example, try:`):
print(`LoadedDie(3,3,{}); `):

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


#Aitken(L): The conjectured limit of the sequence L
#from its three last elements, using Aitken's fromula
Aitken:=proc(L)
local n,a,b,c:
n:=nops(L):

if n<3 then
 RETURN(FAIL):
fi:

a:=L[n-2]: b:=L[n-1]: c:=L[n]:

c-(c-b)^2/(c-2*b+a):
end:


ez:=proc() : print(` Yeladim(v,n), DMP1(m,n), DMP1e(m,n),SeqN(N)  `): 
print(`SeqmN(m,N), SeqmNs(m,N), SP1(m,k), SP(m) `):
print(`WtBlockxq(S,x,q), WtSPxq(S,x,q) `):
print(`GFmxq(m,x,q), DMP1eF(m,n),`):
print(`WtBlockq(S,q) , WtSPq(S,q), GFmq(m,q)`):
print(`SPg(S), GFSq(S,q), DMPS(S,n) `):
print(`MtoVg(M,a,S), SeqSN(S,N)`):
end:


#Yeladim(v,n): given a partition [a1,a2, ..., ak] in
#frequency notation with a1+2*a2+...*k*ak<=n finds
#all the legal continuations to [a1,a2, ..., ak,a_{k+1}] 
#such that no non-zero multiplicity is repeated.
#For example, try:
#Yeladim([1,2],10);
Yeladim:=proc(v,n) local gu,m,i,lu,i1:
m:=add(v[i]*i,i=1..nops(v)):

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

lu:={op(v)} minus {0}:

gu:={[op(v),0]}:

for i1 from 1 to trunc((n-m)/(nops(v)+1)) do 
 if not member(i1,lu) then
  gu:=gu union {[op(v),i1]}:
 fi:
od:

gu:

end:




#DMP1(m,n): the set all distinct-multipliticy partitions with <=n
#with largest part m. For example, try:
#DMP1(3,3);
DMP1:=proc(m,n) local i,mu,mu1:
option remember:

if m=1 then
 RETURN({seq([i],i=0..n)}):
fi:

mu:=DMP1(m-1,n):

{seq(op(Yeladim(mu1,n)),mu1 in mu)}:

end:


SeqNs:=proc(N) local i:
[seq(nops(DMP1(i+1,i+1))-nops(DMP1(i,i)),i=1..N)]:
end:


#SeqmNs(m,N): the first N terms of the sequence
#"partitions of n into parts  <=m with distinct multiplicities"
#slow version
SeqmNs:=proc(m,N) local i:
[1,1,seq(nops(DMP1(m,i+1))-nops(DMP1(m,i)),i=1..N-1)]:
end:

#SeqmN(m,N): the first N terms of the sequence
#"partitions of n into parts  <=m with distinct multiplicities"
SeqmN:=proc(m,N) local i,gu:
gu:=GFmq(m,q):
[seq(coeff(taylor(gu,q=0,N+1),q,i),i=0..N)]:
end:


#SP1(m,k): The set of set partitions of {1, ...,m} into
#k non-empty sets
SP1:=proc(m,k) local gu,mu1,mu,i,gu1:
if m=1 then
 if k=1 then
    RETURN({{{1}}}):
  else
   RETURN({}):
 fi:
fi:

if k=1 then
 RETURN({{{seq(i,i=1..m)}}}):
fi:

gu:=SP1(m-1,k-1):

gu:={seq(gu1 union {{m}},gu1 in gu)}:

mu:=SP1(m-1,k):

for mu1 in mu do

 for i from 1 to nops(mu1) do
  gu:=gu union { {op(1..i-1,mu1),mu1[i] union {m},op(i+1..nops(mu1),mu1)}}:
 od:

od:

gu:

end:


#SP(m): The set of set partitions of {1, ...,m} 
SP:=proc(m) local k:
{seq(op(SP1(m,k)),k=1..m)}:
end:

#WtBlockxq(S,x,q): given a set of integers S finds
#its weight in terms of the indexed variables x and of q
WtBlockxq:=proc(S,x,q) local s:

if nops(S)=1 then
s:=S[1]:
 RETURN(1/(1-q^s*x[s])):
fi:

(-1)^(nops(S)-1)*(nops(S)-1)!*mul(q^s*x[s],s in S)/
(1-mul(q^s*x[s],s in S)):

end:


#WtSPxq(S,x,q): given a set-partition S, finds it weight
#Try:WtSPxq({{1,2}},x,q);
WtSPxq:=proc(S,x,q) local b:
mul(WtBlockxq(b,x,q), b in S)
end:


#DMP1e(m,n): the set of distinct-multiplicity partitions of
#n using {1, ...,m}
DMP1e:=proc(m,n)
DMP1(m,n) minus DMP1(m,n-1):
end:

#GFmxq(m,x,q): The full generating function for integer partitions
#using {1, ...,m} in terms of x[1], ..., x[m] and q
GFmxq:=proc(m,x,q) local gu,gu1:
option remember:
gu:=SP(m):

add(WtSPxq(gu1,x,q),gu1 in gu):

end:

#MtoV(M,a,n): Given a monomial in a[1], ..., a[n]
#finds its list of exponents, for example try
#MtoV(a[1]^3*a[2],3);
MtoV:=proc(M,a,n) local i:
[seq(degree(M,a[i]),i=1..n)]:
end:



#DMP1eF(m,n): the set of distinct-multiplicity partitions of
#n using {1, ...,m} using generating functions
DMP1eF:=proc(m,n) local x,q,gu,i:
gu:=GFmxq(m,x,q):

gu:=taylor(gu,q=0,n+1):

gu:=expand(coeff(gu,q,n)):

{seq(MtoV(op(i,gu),x,m),i=1..nops(gu))}:

end:




#WtBlockq(S,q): given a set of integers S finds
#its weight in terms of q
WtBlockq:=proc(S,q) local s:

if nops(S)=1 then
s:=S[1]:
 RETURN(1/(1-q^s)):
fi:

(-1)^(nops(S)-1)*(nops(S)-1)!*mul(q^s,s in S)/
(1-mul(q^s,s in S)):

end:


#WtSPq(S,q): given a set-partition S, finds it weight
#Try:WtSPq({{1,2}},q);
WtSPq:=proc(S,q) local b:
mul(WtBlockq(b,q), b in S)
end:


#GFmq(m,q): The full generating function for integer partitions
#using {1, ...,m} in terms of q
GFmq:=proc(m,q) local gu,gu1:
option remember:
gu:=SP(m):
add(WtSPq(gu1,q),gu1 in gu):
end:



####start generalized stuff with an arbitrary set S

#SPg(S): The set of set partitions on the list S
SPg:=proc(S) local m,i,lu:
m:=nops(S):
lu:=SP(m):
subs({seq(i=S[i],i=1..m)},lu):
end:


#GFSq(S,q): The full generating function for integer partitions
#with distinct multiplicities
#using the members of S as parts, in terms of q
#For example, try:
#GFSq([1,5,10,25],q);
GFSq:=proc(S,q) local gu,gu1:
option remember:
gu:=SPg(S):

add(WtSPq(gu1,q),gu1 in gu):

end:


#GFSxq(S,x,q): The full generating function for integer partitions
#using the members of S in terms of x[S[1]], ..., x[S[m]] and q
GFSxq:=proc(S,x,q) local gu,gu1:
option remember:
gu:=SPg(S):

add(WtSPxq(gu1,x,q),gu1 in gu):

end:


#MtoVg(M,a,S): Given a monomial M in a[1], ..., a[n]
#finds its list of exponents, [a1,i1] for example try
#MtoVg(a[1]^3*a[2],a,{1,2});
MtoVg:=proc(M,a,S) local i,gu:

gu:=[]:
for i from 1 to nops(S) do
 if degree(M,a[S[i]])>0 then
  gu:=[op(gu),[S[i],degree(M,a[S[i]])]]:
 fi:
od:
gu:
 
end:

#DMPS(S,n): the set of distinct-multiplicity partitions of
#n using the members of S using generating functions
#try DMPS([1,3,4],10);
DMPS:=proc(S,n) local x,q,gu,i:
gu:=GFSxq(S,x,q):

gu:=taylor(gu,q=0,n+1):

gu:=expand(coeff(gu,q,n)):

{seq(MtoVg(op(i,gu),x,S),i=1..nops(gu))}:

end:


#SeqSN(S,N): the first N terms of the sequence
#"partitions of n into parts  drawn from the list S
#with distinct multiplicities"
SeqSN:=proc(S,N) local i,gu,q:
gu:=GFSq(S,q):
[seq(coeff(taylor(gu,q=0,N+1),q,i),i=0..N)]:
end:


###Begin taken from PARTITIONS
#pn(n): the number of partitions of n using Euler's recurrence
#for example, try pn(100);
pn:=proc(n) local j,lu:
option remember:

if n<0 then
RETURN(0):
elif n=0 then
RETURN(1):
else
lu:=0:
for j from 1 while j*(3*j-1)/2<=n do
 lu:=lu-(-1)^j*pn(n-j*(3*j-1)/2):
od:
for j from 1 while j*(3*j+1)/2<=n do
 lu:=lu-(-1)^j*pn(n-j*(3*j+1)/2):
od:
RETURN(lu):
fi:
end:


#pnSeq(n0): the first n0 terms of pn(i). For example, try:
#pnSeq(100);
pnSeq:=proc(n0) local n:
[seq(pn(n),n=1..n0)]:
end:
###End taken from PARTITIONS


#qnmS(n,m,S): the number of partitions of n with largest part <=m
#with distinct multiplicities, none of which belongs to the set S
#For example, try:
#qnmS(10,10,{});
qnmS:=proc(n,m,S) local i,i1:
option remember:
if m=1 then
 if member(n,S) then
   0
 else
   1:
 fi:
else
qnmS(n,m-1,S)+add(qnmS(n-m*i,m-1,S union {i}), 
i in {seq(i1,i1=1..trunc(n/m))} minus S):
fi:
end:

#The number of partitions of n with distinct multiplicity
qn:=proc(n) : qnmS(n,n,{}):end:

#The first N terms of the sequence 
#"number of partitions of n with distinct multiplicity"
qnSeq:=proc(N) local n: [seq(qn(n),n=1..N)]: end:

#The first N terms of the sequence 
#"number of partitions of n with distinct multiplicity"
#using qnPC
qnSeqPC:=proc(N) local n: [seq(qnPC(n),n=1..N)]: end:




 
#SipurDMP(M,q,L): inputs an integer N and outputs
#a webbook listing the (ordinary) generating functions
#for the sequence "number of partitions of n
#with distinct multiplicity using {1, ...,m} for
#m from 1 to M. It also lists the first L terms
#of each sequence
#For example, try:
#SipurDMP(5,q,50);
SipurDMP:=proc(M,q,L) local m,gu,t0,i,n:

t0:=time():
print(`Generating Functions for the Number of Restricted Partitions`):
print(`with Distinct Multiplicities `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):


for m from 1 to M do
gu:=normal(GFmq(m,q)):
print(`Proposition Number`, m):
print(`Let a(n) be the number of partitions of n whose largest part is <=`,m):
print(`and such that all the parts that show up at least once all have`):
print(`different multiplicities. Then`):
print(Sum(a(n)*q^n,n=0..infinity)=gu):
print(`In Maple input format this is:`):
lprint(gu):
print(`For the sake of Sloane, here are the first `, L, `terms of the sequence`):
gu:=taylor(gu,q=0,L+1):
print([seq(coeff(gu,q,i),i=1..L)]):

od:

print(`The webbook was generated in`, time()-t0, `seconds.`):

end:



#SipurChange(S,q,L): inputs a list (or set) of positive integers
# and outputs
#a webbook listing the (ordinary) generating function
#for the sequence "number of partitions of n
#with distinct multiplicity using the elements of S as parts
#It also lists the first L terms
#of each sequence
#For the sake of comparison it also does the (trivial) analogous
#thing when there is no restiction on the multipilicity
#For example, try:
#SipurChange([1,5,10,25],q,50);
SipurChange:=proc(S,q,L) local gu,t0,i,n,Li,Li1,s:

t0:=time():
print(`The Generating Function for the Number of Ways of Making Change`):
print(`of n money-units `):
print(`in a County whose coins have denominatons`, S):
print(`and where all coins (that show up) show up `):
print(` different number of times. `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Prop.: Let a(n) be the number of ways of forming change`):
print(`of n money-units in a country where the coins have denominations`):
print(S):
print(`and none of the coins that do show up at least once `):
print(` can occur the same number of`):
print(`times. `):
gu:=GFSq(S,q):
print(`Then `):
print(Sum(a(n)*q^n,n=0..infinity)=gu):
print(`In Maple input format this is:`):
lprint(gu):
gu:=taylor(gu,q=0,L+1):
Li:=[seq(coeff(gu,q,i),i=1..L)]:

print(`For the sake of Sloane, the first`, L, `terms are`):
print(Li):

gu:=mul(1/(1-q^s), s in S):

print(`For the sake of comparison the generating function for`):
print(`the number of ways of making change without restiction is,`):
print(`of course, the much simpler`):
print(gu):

gu:=taylor(gu,q=0,L+1):

Li1:=[seq(coeff(gu,q,i),i=1..L)]:

print(`For the sake of Sloane, the first`, L, `terms are`):
print(Li1):

print(`The sequence of ratios of good/all is`):
print(evalf([seq(Li[i]/Li1[i],i=1..L)])):

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

print(`This took`, time()-t0, `seconds .`):

end:




#SipurSeq(L): inputs an integer L 
# and outputs the first L terms of the
#sequence "number of partitions of n with distinct multiplicity"
#It also lists the first L terms
#of the classical sequence "number of partitions"
#and tries to find the asymptotic trend
#SipurSeq(100);
SipurSeq:=proc(L) local gu,gu1,i,t0:

t0:=time():
gu:=qnSeq(L):
print(`The First `, L, `terms of the sequence the number of`):
print(` OEIS sequence A098859  `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`The first L terms of the sequence:`):
print(` "number of distinct multiplicity partitions of n" are:`):
print(gu):
print(`For the sake of comparison, the sequence for the number of good all`):
print(`partitions is`):
gu1:=pnSeq(L):
print(gu1):
print(`The sequence of ratios good/all is:`):
print(evalf([seq(gu[i]/gu1[i],i=1..L)])):
print(`Going back to the first, interesting, sequence`):
print(`the sequence of log(a(n))/sqrt(n) is`):
print(evalf([seq(log(gu[i])/sqrt(i),i=1..L)])):
print(`that gives some support to the conjecture that the asymptotics`):
print(`is roughly exp(C*sqrt(n)) for some constant C, that we have no`):
print(`idea what it is, but it should be less than`):
print(Pi*sqrt(2/3)=evalf(Pi*sqrt(2/3))):
print(`that, according to Ramanujan is the constant`):
print(`for the classical paritions.`):
print(``):
print(`The whole thing took`, time()-t0, `seconds. `):

end:



#SipurSeqPC(L): inputs an integer L 
# and outputs the first L terms of the
#sequence "number of partitions of n with distinct multiplicity"
#using precomputed values
#It also lists the first L terms
#of the classical sequence "number of partitions"
#and tries to find the asymptotic trend
#SipurSeqPC(100);
SipurSeqPC:=proc(L) local gu,gu1,i,t0:

if L>508 then
print(`L must be <=508`):
RETURN(FAIL):
fi:

t0:=time():
gu:=qnSeqPC(L):
print(`The First `, L, `terms of the sequence the number of`):
print(` OEIS sequence A098859  `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`The first L terms of the sequence:`):
print(` "number of distinct multiplicity partitions of n" are:`):
print(gu):
print(`For the sake of comparison, the sequence for the number of good all`):
print(`partitions is`):
gu1:=pnSeq(L):
print(gu1):
print(`The sequence of ratios good/all is:`):
print(evalf([seq(gu[i]/gu1[i],i=1..L)])):

print(`the sequence of log(a(n))/sqrt(n) for the partition function is`):
print(evalf([seq(log(gu1[i])/sqrt(i),i=1..L)])):

print(`We know thanks to Hardy and  Ramanujan that this converges to`):
print(Pi*sqrt(2/3)=evalf(Pi*sqrt(2/3))):
print(`--------------------------------------------------`):
print(`Going back to the first, interesting, sequence`):
print(`the sequence of log(a(n))/sqrt(n) is`):
print(evalf([seq(log(gu[i])/sqrt(i),i=1..L)])):

print(`that gives some support to the conjecture that the asymptotics`):
print(`is roughly exp(C*sqrt(n)) for some constant C, that we have no`):
print(`idea what it is, but it should be less than`):
print(Pi*sqrt(2/3)=evalf(Pi*sqrt(2/3))):
print(`that, according to Ramanujan is the constant`):
print(`for the classical paritions.`):
print(``):
print(`The whole thing took`, time()-t0, `seconds. `):

end:





#qnPC(n): The precomputed value of qn(n) using Maciej Ireneusz Wilczynsk's
#table from the OEIS for the first 508 values
qnPC:=proc(n) local Seq508:
Seq508:=
[ 1,
 2,
 2,
 4,
 5,
 7,
 10,
 13,
 15,
 21,
 28,
 31,
 45,
 55,
 62,
 82,
 105,
 116,
 153,
 172,
 208,
 251,
 312,
 341,
 431,
 492,
 588,
 676,
 826,
 905,
 1120,
 1249,
 1475,
 1676,
 2003,
 2187,
 2625,
 2922,
 3409,
 3810,
 4481,
 4910,
 5792,
 6382,
 7407,
 8186,
 9527,
 10434,
 12104,
 13308,
 15293,
 16825,
 19332,
 21193,
 24281,
 26590,
 30264,
 33167,
 37746,
 41304,
 46723,
 51105,
 57824,
 63089,
 71116,
 77661,
 87174,
 95162,
 106767,
 116276,
 130086,
 141975,
 158109,
 172158,
 192196,
 208750,
 231985,
 252703,
 279917,
 304481,
 337512,
 365938,
 404723,
 440290,
 484899,
 526012,
 580656,
 628704,
 691427,
 750854,
 823135,
 891834,
 979467,
 1058435,
 1159767,
 1256880,
 1371965,
 1483751,
 1623375,
 1752443,
 1911046,
 2067456,
 2249444,
 2429337,
 2647532,
 2852449,
 3101167,
 3350292,
 3632299,
 3916575,
 4254015,
 4578913,
 4961281,
 5350317,
 5786186,
 6229771,
 6744992,
 7249767,
 7834973,
 8438264,
 9099153,
 9782250,
10564186,
11342088,
 12225064,
 13143216,
 14142579,
 15185699,
 16355999,
 17537617,
 18859435,
 20249551,
 21740304,
 23309862,
 25052839,
 26830676,
 28789722,
 30864369,
 33075342,
 35421402,
 37984022,
 40629214,
 43513946,
 46590045,
 49835545,
 53295668,
 57045908,
 60952531,
 65159723,
 69668751,
 74401163,
 79476261,
 84912875,
 90619911,
 96721671,
 103287973,
 110130121,
 117497296,
 125334380,
 133622042,
 142393866,
 151868166,
 161705680,
 172335955,
 183548697,
 195467141,
 208011607,
 221601927,
 235641300,
 250847732,
 266797381,
 283843814,
 301654570,
 320989996,
 340913983,
 362541686,
 385082756,
 409262515,
 434418643,
 461778184,
 489866305,
 520392683,
 552079521,
 586194919,
 621502624,
 659937596,
 699321044,
 742181954,
 786463956,
 834255416,
 883555471,
 937267251,
 992173017,
 1051948501,
 1113519332,
 1180126172,
 1248576566,
 1323165871,
 1399318823,
 1482270957,
 1567435521,
 1659678555,
 1754248355,
 1857336621,
 1962414004,
 2076850407,
 2194072607,
 2321208607,
 2451231508,
 2592891255,
 2737192203,
 2894320440,
 3054910670,
 3229181482,
 3407102533,
 3600920587,
 3798138044,
 4012753123,
 4231776835,
 4469573493,
 4711985790,
 4975800649,
 5244141186,
 5536011355,
 5833446461,
 6156363497,
 6485152020,
 6842757583,
 7206338990,
 7601390509,
 8003621810,
 8440255173,
 8884410449,
 9366957023,
 9857479220,
 10390079693,
 10931848047,
 11519618594,
 12117163210,
 12765729427,
 13425002589,
 14139888851,
 14866704845,
 15654865741,
 16455709018,
 17323894915,
 18206361345,
 19162301209,
 20133836524,
 21186368654,
 22255630035,
 23413391563,
 24590375112,
 25863657538,
 27157363736,
 28557916275,
 29980415516,
 31518521838,
 33082462559,
 34772331361,
 36489309069,
 38345824594,
 40231419904,
 42267733469,
 44338792122,
 46573550848,
 48844035616,
 51296555759,
 53787805828,
 56474171182,
 59207337566,
 62152596011,
 65145449003,
 68374138005,
 71654266841,
 75186817149,
 78782126087,
 82650818863,
 86582791946,
 90819450417,
 95125166564,
 99755056896,
 104469008712,
 109534680015,
 114684788687,
 120226625324,
 125860184790,
 131909730528,
 138072023455,
 144684256923,
 151409122264,
 158635914685,
 165985533368,
 173866081523,
 181896839384,
 190502591918,
 199258459499,
 208654515323,
 218214161780,
 228450132468,
 238886873736,
 250054289424,
 261421189189,
 273602929772,
 286002865649,
 299261014944,
 312785051774,
 327236350732,
 341953298978,
 357702466118,
 373741615603,
 390866962134,
 408345191419,
 426995222970,
 445996749126,
 466303861629,
 486995844122,
 509057737739,
 531584882279,
 555590087871,
 580061457660,
 606176667244,
 632800307198,
 661148541551,
 690109829232,
 720929273785,
 752362421550,
 785861882202,
 820032808772,
 856367089563,
 893505180557,
 932974327023,
 973251787878,
 1016117420564,
 1059864741504,
 1106321649251,
 1153830594036,
 1204254772572,
 1255737778694,
 1310456344477,
 1366332301606,
 1425590792679,
 1486219133035,
 1550487103296,
 1616140458601,
 1685826947551,
 1757023692803,
 1832436289827,
 1909630294234,
 1991354959336,
 2074883670289,
 2163429792027,
 2253944402366,
 2349698790863,
 2447758614720,
 2551449289557,
 2657485099480,
 2769745998371,
 2884559965889,
 3005875826937,
 3130170084566,
 3261444850084,
 3395753881732,
 3537772684228,
 3683097037902,
 3836469501420,
 3993673806819,
 4159516599395,
 4329278324570,
 4508565901992,
 4692113934529,
 4885608683026,
 5084022842150,
 5293103607420,
 5507224316064,
 5733094553259,
 5964446261482,
 6208063062312,
 6457971050222,
 6721029131609,
 6990555023173,
 7274543535740,
 7565554524028,
 7871665923014,
 8185807128084,
 8516125822850,
 8854715520640,
 9211076581173,
 9576414686942,
 9960308958046,
 10354413887773,
 10768393301685,
 11192922947428,
 11639248496782,
 12097009929995,
 12577535626851,
 13071023045341,
 13588878870627,
 14120145976016,
 14678104663727,
 15250594069993,
 15850968802464,
 16467739212256,
 17114343867415,
 17777960310607,
 18474202653142,
 19188856253478,
 19937610272446,
 20707061191920,
 21512980271094,
 22340403967845,
 23207658190803,
 24098177962752,
 25030335001617,
 25988544761734,
 26991259939139,
 28021108499306,
 29099499600709,
 30207205055458,
 31365679647039,
 32556869122455,
 33802307171692,
 35081871490402,
 36420511224176,
 37796006046691,
 39233306527839,
 40711597424928,
 42255905533271,
 43843038722757,
 45501973060768,
 47207101834491,
 48987388949023,
 50818908129487,
 52730662431804,
 54696022041129,
 56748512690792,
 58858799744551,
 61060344050382,
 63325771382121,
 65688586609411,
 68118330144078,
 70653706224513,
 73261190414204,
 75979336025739,
 78776976025539,
 81692658854391,
 84691739975535,
 87818687188342,
 91035402516895,
 94386160017132,
 97835594272348,
 101427958179676,
 105123981932355,
 108974616938236,
 112936758662400,
 117061022106221,
 121307567966840,
 125726951995918,
 130274899715040,
 135009572015684,
 139882418241607,
 144951192623635,
 150171149219672,
 155599852452603,
 161187678160663,
 167000737395708,
 172984720826185,
 179205147001419,
 185612115322290,
 192270957480172,
 199126292431341,
 206253058424261,
 213590756226113,
 221213517690183,
 229066033690302,
 237222036228622,
 245620259808715,
 254345213656589,
 263329980201515,
 272658034409548,
 282268446864145,
 292244288373406,
 302518080164827,
 313184727002448,
 324170760655573,
 335569752459209,
 347315247528261,
 359500163053252,
 372050860369759,
 385073528287008,
 398488027340602,
 412398846972096,
 426733987212084,
 441597072729282,
 456908448135995,
 472786123885576,
 489143704939202,
 506097184174830,
 523569441563756,
 541675325667328,
 560329445596724,
 579662655673247,
 599582531634653,
 620217094197502,
 641484786562528,
 663512198980723,
 686209157979683,
 709719461956859]:

if not type(n,integer) or n<1 then
print(`Bad input`):
RETURN(FAIL):
fi:

if n>508 then
print(`pre-computed values are only available up to 508`):
RETURN(FAIL):
fi:

Seq508[n]:

end:






#LoadedDie(n,m,S): picks the multiplicity of m in a random
#distinct-multiplicity partition of n with largest part
#<=m and avoiding multiplicities in S. For example, try:
#LoadedDie(n,m,S) 
LoadedDie:=proc(n,m,S) local lu,i:

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

lu:=[]:
if qnmS(n,m-1,S)>0 then
 lu:=[op(lu),[[n,m-1,S,0],qnmS(n,m-1,S)]]:
fi:

for i from 1 to trunc(n/m) do
 if not member(i,S) and qnmS(n-m*i,m-1,S union {i})>0 then
   lu:=[op(lu),[[n-m*i,m-1,S union {i},i],qnmS(n-m*i,m-1,S union {i})]]:
 fi:
od:


if add(lu[i][2],i=1..nops(lu))<>qnmS(n,m,S) then
 ERROR(`something is wrong`):
fi:

lu:

end:





#RandDMP1(n,m,S): A uniformly-at-random distinct-multiplicity
#partition of n with largest part <=m
#none of which belongs to the set S
#For example, try:
#RandDMP1(10,10,{});
RandDMP1:=proc(n,m,S) local lu,N,r,j,P,i,mu,lu1:

if n=0 then
 RETURN([]):
fi:

if m=1 then
 if not member(n,S) then
    RETURN([[1,n]]):
  else
   RETURN(FAIL):
 fi:
fi:

lu:=LoadedDie(n,m,S):

mu:=[seq(lu[i][2],i=1..nops(lu))]:

N:=convert(mu,`+`):

r:=rand(1..N)():

for i from 1 to nops(mu) while add(mu[j],j=1..i) <r do od:


lu1:=lu[i][1]:

P:=RandDMP1(op(1..3,lu1)):
if lu1[4]=0 then
 RETURN(P):
else
 RETURN([op(P),[m,lu1[4]] ]):
fi:
end:




#RandDMP(n): A uniformly-at-random distinct-multiplicity
#partition of n 
#For example, try:
#RandDMP(10);
RandDMP:=proc(n) RandDMP1(n,n,{}): end: