######################################################################
##FCP.txt: Save this file as  PCP.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read PCP.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Oct. 17, 2018
#This version: Nov. 6, 2018 (to add procedure Ryba(R,x) implementing Chris Ryba's brilliant solution of Question 2
print(`Created:  Oct. 17, 2018`):
print(` This is PCP.txt `):
print(`It is a Maple package that accompanies the article `):
print(`Fractional Counting of Integer Partitions`):
print(`by  Doron Zeilberger and Noam Zeilberger`):

print(`available from Doron Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.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(plots):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: bnBF, bnkGF, Fr, PlotRyba, SeqbnGF `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  bn, bnF, bnk, bnkF, NoamC, PlotB, Ryba, Seqbn`):
 print(` `):

elif nops([args])=1 and op(1,[args])=bn then
print(`bn(n): the fractional count of theset of integer partitions of n done by dynamical programming, hence much faster. `):
print(`Try: bn(10);`):

elif nops([args])=1 and op(1,[args])=bnF then
print(`bnF(n): the fractional count of theset of integer partitions of n done by dynamical programming in Floating point, hence much faster. `):
print(`Try: bnF(10);`):

elif nops([args])=1 and op(1,[args])=bnk then
print(`bnk(n,k): the sum of 1/productOfParts over all partitions of n whose largest part is k. Try:`):
print(`bnk(30,20);`):

elif nops([args])=1 and op(1,[args])=bnkF then
print(`bnkF(n,k): like bnk(n,k) (q.v.) but in floating point`):
print(`Try:`):
print(`bnkF(30,20);`):

elif nops([args])=1 and op(1,[args])=bnkGF then
print(`bnkGF(n,k): like bnk(n,k) (q.v.) but via generating function. Try: `):
print(`bnkGF(30,20);`):

elif nops([args])=1 and op(1,[args])=bnBF then
print(`bnBF(n): the fractional count of theset of integer partitions of n done by brute force. Very slow. Only for checking purposes. `):
print(`should be the same as bn(n); `):
print(`Try: bnBF(10);`):

elif nops([args])=1 and op(1,[args])=Fr then
print(`Fr(r,x): the function F_r(x) defined by Chris Ryba (see his paper). Try:`):
print(`Fr(2,x);`):

elif nops([args])=1 and op(1,[args])=NoamC then
print(`NoamC(): the Last 11 precomputed values of  SeqbnN(15000)`):

elif nops([args])=1 and op(1,[args])=PlotB then
print(`PlotB(n): The plot of [seq([k/n,bnkF(n,k)],k=1..n)]). Try:`):
print(`PlotB(1000);`):

elif nops([args])=1 and op(1,[args])=PlotRyba then
print(`PlotRyba(R,x,N): Plots the first R members of Chris Ryba's functions F_r(x). It also plots bnkF(N,x*N) for that range;`):
print(`R must be 1,2, or 3. Try: `):
print(`PlotRyba(3,x,1000);`):

elif nops([args])=1 and op(1,[args])=Ryba then
print(`Ryba(R,x): the first R members of Christopher Ryba's functions F_r(x). Try:`):
print(`Ryba(4,x);`):

elif nops([args])=1 and op(1,[args])=Seqbn then
print(`Seqbn(N): the sequence of bn(n)  for n from 1 to N. Try:`):
print(`Seqbn(100);`):

elif nops([args])=1 and op(1,[args])=SeqbnGF then
print(`SeqbnGF(N): Like Seqbn(N) but via the generating function. Try: `):
print(`SeqbnGF(10);`):

elif nops([args])=1 and op(1,[args])=SeqbnN then
print(` SeqbnN(N): the sequence of bn(n)/n for n from 1 to N. Do:`):
print(`SeqbnN(1000);`):

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



ez:=proc(): print(`  EstC(L,n,k), EstCg(L,n,k,alpha), Aitken(L)  `): 
print(`Gkz(3,z), CheckGkz(k,N), PlotB(n) , EstNC(n,k,r) `):
end:

with(combinat):

#bnBF(n): the fractional count of theset of integer partitions of n done by brute force. Try: bnBF(10);
bnBF:=proc(n) local gu,gu1,i1:
gu:=partition(n):
add(1/mul(gu1[i1],i1=1..nops(gu1)),gu1 in gu):
end:

#bnk(n,k): the sum of 1/productOfParts over all partitions of n whose largest part is k
bnk:=proc(n,k) option remember:
if k=1 then
 RETURN(1):
fi:

if k=n then
 RETURN(1/n):
fi:

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

(k-1)/k*bnk(n-1,k-1)+1/k*bnk(n-k,k):
end:

bn:=proc(n) local k:
add(bnk(n,k),k=1..n):
end:

#Seqbn(N): the sequence of bn(n)  for n from 1 to N
Seqbn:=proc(N) local n:
[seq(bn(n),n=1..N)]:
end:


#SeqbnN(N): the sequence of bn(n)/n
SeqbnN:=proc(N) local n:
[seq(bn(n)/n,n=1..N)]:
end:



EstCg:=proc(L,n,k,alpha) local gu,a,n1,i,var,eq:
if not k<nops(L)/10 then
  RETURN(FAIL):
fi:

gu:=add(a[i]/n^(alpha*i),i=0..k):
var:={seq(a[i],i=0..k)}:
eq:={seq(subs(n=n1,gu)-L[n1],n1=nops(L)-k..nops(L))}:
var:=solve(eq,var):
subs(var,gu):
end:

##########Float



#bnkF(n,k): the sum of 1/productOfParts over all partitions of n whose largest part is k
bnkF:=proc(n,k) option remember:
if k=1 then
 RETURN(1):
fi:

if k=n then
 RETURN(evalf(1/n)):
fi:

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

evalf((k-1)/k*bnkF(n-1,k-1)+1/k*bnkF(n-k,k)):
end:

bnF:=proc(n) local k:
add(bnkF(n,k),k=1..n):
end:

#SeqbnF(N): the sequence of b(n)  for n from 1 to N
SeqbnF:=proc(N) local n:
[seq(bnF(n),n=1..N)]:
end:


#SeqbN(N): the sequence of b(n)/n
SeqbnFN:=proc(N) local n:
[seq(bnF(n)/n,n=1..N)]:
end:

#Aitken(Lis1): The Aitken transform of the sequence Lis1
Aitken:=proc(Lis1)
local gu,n:
gu:=[]:

for n from 3 to nops(Lis1) do
gu:=[op(gu),(Lis1[n]*Lis1[n-2]-Lis1[n-1]^2)/(Lis1[n]-2*Lis1[n-1]+Lis1[n-2])]:
od:

gu:
end:

#Gkz(k,z): The generating function Sum(bnk(n,k)*z^n,n=0..infinity); Try:
#Gkz(3,z):
Gkz:=proc(k,z) local i:

z^k/k/mul(1-z^i/i, i=1..k):

end:


#CheckGkz(k,N) local gu,z:
CheckGkz:=proc(k,N) local gu,z,i:
gu:=Gkz(k,z):
gu:=taylor(gu,z=0,N+1):
evalb({seq(coeff(gu,z,i)-bnk(i,k),i=0..N)}={0}):
end:

#PlotB(n): The plot of [seq([k/n,bnkF(n,k)],k=1..n)]). Try:
#PlotB(1000);
PlotB:=proc(n) local k:
plot([seq([k/n,bnkF(n,k)],k=1..n)]):
end:


#NoamC(): the Last 11 precomputed values of  SeqbnN(15000)
NoamC:=proc()
[0.5611411658, 0.5611411846, 0.5611412033, 0.5611412220, 0.5611412407,

    0.5611412594, 0.5611412781, 0.5611412968, 0.5611413156, 0.5611413344,

    0.5611413530]:
end:



#EstC(L,n,k,N): given a list L whose last entry is the value of a sequence at N tries to guess
#the asymptotic
EstC:=proc(L,n,k,N) local gu,a,n1,i,var,eq:

gu:=add(a[i]/n^i,i=0..k):
var:={seq(a[i],i=0..k)}:
eq:={seq(subs(n=N-n1,gu)-L[nops(L)-n1],n1=0..k)}:
var:=solve(eq,var):
subs(var,gu):
end:

EstNC:=proc(n,k,r):
if r>=12 or r<=k+1 then
  RETURN(FAIL):
fi:
EstC([op(1..r,NoamC())],n,k,15000-11+r):
end:



#bnkFg(n,k,alpha): the sum of 1/(productOfParts)^alpha over all partitions of n whose largest part is k
bnkFg:=proc(n,k,alpha) option remember:
if k=1 then
 RETURN(1):
fi:

if k=n then
 RETURN(evalf(1/n^alpha)):
fi:

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

evalf(((k-1)/k)^alpha*bnkFg(n-1,k-1,alpha)+1/k^alpha*bnkFg(n-k,k,alpha)):
end:

bnF:=proc(n) local k:
add(bnkF(n,k),k=1..n):
end:

bnFg:=proc(n,alpha) local k:
add(bnkFg(n,k,alpha),k=1..n):
end:

#SeqbnFg(N,alpha): the sequence of b(n)  for n from 1 to N
SeqbnF:=proc(N,alpha) local n:
[seq(bnFg(n,alpha),n=1..N)]:
end:

#PlotBg(n,alpha): The plot of [seq([k/n,bnkF(n,k)],k=1..n)]). Try:
#PlotBg(1000,alpha);
PlotBg:=proc(n,alpha) local k:
plot([seq([k/n,bnkFg(n,k,alpha)],k=1..n)]):
end:

#SeqbnGF(N): Like Seqbn(N) but via the generating function
SeqbnGF:=proc(N) local gu,i,q:
gu:=mul(1/(1-q^i/i),i=1..N):
gu:=taylor(gu,q=0,N+1):
[seq(coeff(gu,q,i),i=1..N)]:
end:


#bnkGF(n,k): Like bnk(n,k) but via generating functions
bnkGF:=proc(n,k) local gu,i,q:
gu:=q^k/k*mul(1/(1-q^i/i),i=1..k):
gu:=taylor(gu,q=0,n+1):
coeff(gu,q,n):
end:


MyInt:=proc(f,x,a,b) local gu:
gu:=int(f,x):
subs(x=b,gu)-subs(x=a,gu):
end:

#Fr(r,x): the function F_r(x) defined by Chris Ryba (see his paper). Try:
#Fr(2,x);
Fr:=proc(r,x) local s,t:
option remember:
if (not type(r,integer) and r>=1 and type(x,symbol)) then
print(`bad input`):
RETURN(FAIL):
fi:
if r=1 then
 RETURN((1-x)/x):
fi:

(1-x)/x-(1-x)/x*
(
MyInt(subs(x=t,Fr(r-1,x)),t,x/(1-x),1/(r-1))
+add(
MyInt(subs(x=t,Fr(s,x)),t,1/(s+1),1/s),
s=1..r-2)
)
:

end:


#Ryba(R,x): the first R members of Chris Ryba's functions F_r(x). Try:
#Ryba(4,x);
Ryba:=proc(R,x) local r:


[seq(Fr(r,x),r=1..R)]:
end:

#PlotRyba(R,x,N): Plots the first R members of Chris Ryba's functions F_r(x). Try:
#PlotRyba(4,x,N);
PlotRyba:=proc(R,x,N) local k,r,gu,pic,pic1:
if not (type(R, integer) and R>=1 and R<=3) then
 print(`R must be 1, 2, or 3`):
 RETURN(FAIL):
fi:

gu:=Ryba(R,x):
gu:=evalf(exp(-gamma)*gu):
pic:=plot(gu[1],x=1/2..1,color=red):

for r from 2 to R do
pic:=pic,plot(gu[1],x=1/(r+1)..1/r,color=red):
od:

pic1:=plot([seq([k/N,bnkF(N,k)],k=trunc(N/(R+1))..N)],color=blue):

pic:=pic,pic1:
display(pic);

end: