######################################################################
## DymLuks.txt Save this file as   DymLuks.txt to use it,            #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `DymLuks.txt` <Enter>                        #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################



print(`First Written: Jan. 2023: tested for Maple 2020 `):
print(`This Version  : Jan. 2, 2023  `):
print():
print(`This is  DymLuks.txt, a Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` Experimenting with the Dym-Luks Ball and Cell Game (almost) Sixty Years Later `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/DymLuks.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the Story functions type: ezraSt();`):
print(`For a list of the SIMULTATION functions type: ezraSi();`):
print(`For a list of the Checking functions type: ezraC();`):
print(`For a list of the General functions type: ezraG();`):
print():





ezraG:=proc()
if args=NULL then

print(`The General procedures are: GFg, GuessExpPol1 ,GuessExpoPol, STATg, SymMoms, SymMomsS, SymMomsSlim  `):

else
ezra(args):

fi:

end:


ezraSi:=proc()
if args=NULL then

print(`The Simulation procedures are: DL, DLv, OS, OSv, RF, SimuGFrn, SimuSTATrn `):

else
ezra(args):

fi:

end:

ezraSt:=proc()
if args=NULL then

print(`The Story procedures are: Paper1, Paper1t, Paper2, Paper2t, Paper3, Paper3t, Paper4, Paper4R, Paper5 `):

else
ezra(args):

fi:

end:

ezraC:=proc()
if args=NULL then

print(`The Checking procedures are:   FUNrn, TPbf `):

else
ezra(args):

fi:


end:


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: EstC, GFrnA, Mnr, M2nr, NuCap, Prt, TP, Stat, STATrnA  `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` DymLuks.txt: A Maple package for  Experimenting with the Dym-Luks Ball and Cell game`):

print(`The MAIN procedures are:  Enr, GFrn, LimDL, PlotDist, STATrn  `):
print(``):


elif nargs=1 and args[1]=DL then
print(`DL(r,n): Playing the Dym-Luks Ball and Cell game and recording its duration. Try:`):
print(`DL(10,10);`):


elif nargs=1 and args[1]=DLv then
print(`DLv(r,n): A verbose version of DL(r,n) (q.v.). Try:`):
print(`DLv(10,10);`):

elif nargs=1 and args[1]=Enr then
print(`Enr(n,r): Mnr(n,r)-Sum((1/j)*(n/(n-1))^(j-1),j=1..r): Try: `):
print(`Enr(2,10);`):


elif nargs=1 and args[1]=EnrR then
print(`EnrR(n,r): Mnr(n,r)/Sum((1/j)*(n/(n-1))^(j-1),j=1..r): Try: `):
print(`EnrR(2,10);`):


elif nargs=1 and args[1]=E2nrR then
print(`E2nrR(n,r): M2nr(n,r)/(Sum((1/j)*(n/(n-1))^(j-1),j=1..r)-Sum(((1/j)*(n/(n-1))^(j-1))^2,j=1..r))). Try: `):
print(`E2nr(3,40);`):


elif nargs=1 and args[1]=EstC then
print(`EstC(L,d,r): Given a sequence of numbers L believed to have the asympotic expansion C[0]+C[1]/sqrt(r)+C[2]/r+.. estimates assuming that`):
print(`It returns the list [C[0],.., C[d]], looking at the last d+1 terms. It also returns looking at d+1 terms before that.`):
print(`it ends in d. Try:`):
print(`EstC([seq(5+1/sqrt(i)+3/i,i=1..100)],2,r);`):

elif nargs=1 and args[1]=FUNrn then
print(`FUNrn(r,n): all the functions from {1,..,r} to {1,..,n} given as lists of length r with entries in {1,...,n}. Try:`):
print(`FUNrn(2,2);`):

elif nargs=1 and args[1]=GFg then
print(`GFg(n,a, x,n1): The duration probability generating function for the Markov process with prob. of n->n-1 given by the expression, a(n), and the prob. that n->n  by (1-a(n)). If it is currently at n1. Try:`):
print(`GFg(n,1/2^n,x,10);`):

elif nargs=1 and args[1]=GFrn then
print(`GFrn(r,n,x): the rational function in x giving the probability generating function for the duration of the Dym-Luks Ball and Cell Game with r balls and n cells. Try:`):
print(`GFrn(4,4,x);`):

elif nargs=1 and args[1]=GFrnA then
print(`GFrnA(r,n,x): the rational function that is the approximation to GFrn(r,n,x) for small n and large r. Easier to deal with. Try:`):
print(`GFrnA(4,4,x);`):


elif nargs=1 and args[1]=GuessExpPol then
print(`GuessExpPol(L,d1,d2,a,n): Given a list of expressions in a and a^n, L, finds a polyomial of degree<= d1 in a^n and degree <=d2 in n, let's call it POL such that`):
print(`L[n1]=POL(n1,a^n1). Try:`):
print(`GuessExpPol([seq(n1*a^n1,n1=1..10)],1,1,a,n);`):

elif nargs=1 and args[1]=GuessExpPol1 then
print(`GuessExpPol1(L,d1,d2,a,n): Given a list of expressions in a and a^n, L, finds a polyomial of degree d1 in a^n and degree d2 in n, let's call it POL such that`):
print(`L[n1]=POL(n1,a^n1). Try:`):
print(`GuessExpPol1([seq(n1*a^n1,n1=1..10)],1,1,a,n);`):

elif nargs=1 and args[1]=LimDL then
print(`LimDL(n,R,K): For a fixed (numeric) n, and positive integer R, and positive integer K,`):
print(` finds K lists. The first list consists of the first R terms M_n(r) (the average duration of the Dym-Luks game with r balls and n cells) for r from 2 to R. `):
print(`The second list of the ratio of s.d/average, and the 3rd through K lists the scaled moments. Try:`):
print(`LimDL(3,20,4); `):

elif nargs=1 and args[1]=Mnr then
print(`Mnr(n,r): STATrn(r,n,2)[1], i.e. the average duration of a Dym-Luks game with  r balls and n cells. Try:`):
print(`Mnr(2,10);`):


elif nargs=1 and args[1]=M2nr then
print(`M2nr(n,r): STATrn(r,n,2)[2]^2, i.e. the variance of duration of a Dym-Luks game with  r balls and n cells. Try:`):
print(`M2nr(2,10);`):

elif nargs=1 and args[1]=NuCap then
print(`NuCap(L): Given a list of integers L outputs the number of entries that show up exactly once. Try:`):
print(`nuCap([5,4,2,2,1,1,3]);`):

elif nargs=1 and args[1]=OS then
print(`OS(r,n): One step of the Dym/Luks game if there are r balls, and n cells. It returns the number of balls not captured. Try:`):
print(`OS(10,10);`):

elif nargs=1 and args[1]=OSv then
print(`OSv(r,n): Verbose version of OS(r,n) (q.v.). Try:`):
print(`OSv(100,100);`):

elif nargs=1 and args[1]=Paper1 then
print(`Paper1(R,n,K): A paper about the Dym-Luks Ball and Cell Game with fixed r and general (symbolic) n. For r=1,...R, it gives the`):
print(`prob. generating function and all the moments up to the K's. Try:`):
print(`Paper1(4,n,6);`):

elif nargs=1 and args[1]=Paper1t then
print(`Paper1t(R,n,K): An abbreviated version of Paper1(R,n,K) (q.v.) only the average, s.d. and scaled moments. try:`):
print(`Paper1t(4,n,6);`):


elif nargs=1 and args[1]=Paper2 then
print(`Paper2(R,N,K): A paper about the Dym-Luks Ball and Cell Game with fixed (small) n from 1 to N,  and r from 1 to R, it gives the`):
print(`prob. generating function and all the moments up to the K's. Try:`):
print(`Paper2(5,5,4);`):


elif nargs=1 and args[1]=Paper2t then
print(`Paper2t(R,N,K): Like  Paper2(R,N,K) but only with moments. Try: `):
print(`Paper2t(5,5,4);`):

elif nargs=1 and args[1]=Paper3 then
print(`Paper3(R,a,K): A paper about the Dym-Luks Ball and Cell Game with r balls and trunc(a*r) cells`):
print(`prob. generating function and all the moments up to the K's. Try:`):
print(`Paper3(10,1,4);`):

elif nargs=1 and args[1]=Paper3t then
print(`Paper3t(R,a,K): Like Paper3(R,a,K) (q.v.) but only moments (i.e. no prob. gen. functions). Try:`):
print(`Paper3t(10,1,4);`):

elif nargs=1 and args[1]=Paper4 then
print(`Paper4(N,R): A table of the quantity Enr(n,r) for n from 2 to N and r from 1 to R. Try:`):
print(`Paper4(3,40);`):


elif nargs=1 and args[1]=Paper4R then
print(`Paper4R(N,R): A table of the quantities EnrR(n,r) and E2nrR(n,r), for n from 2 to N and r from 1 to R. Try:`):
print(`Paper4(3,40);`):

elif nargs=1 and args[1]=Paper5 then
print(`Paper5(K): an article about the expectation, variance, and moments up to the K-th (including the limit as n goes to infinity) for the duration`):
print(`of a Markov process from n to 0 where the prob. of moving from n to n-1 is a^n, and of staying at n is 1-a^n. Try:`):
print(`Paper5(4);`):

elif nargs=1 and args[1]=PlotDist then
print(`PlotDist(r,n,K): plots the distribution of the duration of the Dym-Luks Ball and Cell game with r balls and n cells, using up to the K-th Taylor expansion. Try:`):
print(`PlotDist(10,10,30);`):

elif nargs=1 and args[1]=Prt then
print(`Prt(n,r,t): The probability of getting from state r to state r-t, what the Dym-Luks paper calls P_{r,r-t}(n). Try:`):
print(` Prt(2,r,1);`):

elif nargs=1 and args[1]=RF then
print(`RF(r,n): a random placing of balls {1,..,r} into cells {1,...,n} it outputs a list of length r where L[i] is the cell where ball is placed. Try:`):
print(`RF(5,7);`):

elif nargs=1 and args[1]=SimuGFrn then
print(`SimuGFrn(r,n,x,K): an approximation to GFrn(r,n,x) (q.v.) by playing the game K times. Try:`):
print(`Of course, every time you get a (slightly) different output (since this is random). As K gets bigger, of course, `):
print(`you should get more agreement. Try:`):
print(`SimuGFrn(3,3,x,100);`):


elif nargs=1 and args[1]=SimuSTATrn then
print(`SimuSTATrn(r,n,K1,K2): Simulates the Dym-Luks game K1 time and finds the average, s.d., and the third, ..., up to the K2-th moment. Try:`):
print(`SimuSTATrn(10,10,100,6);`):

elif nargs=1 and args[1]=Stat then
print(`Stat(F,t,K): Given a generating function F of t, for a r.v. , and a positive ineger K>=2, finds the statistical information, i.e. the list of length K`):
print(`consisiting of the  expectation, followed by the standard-deviation, forllowed by the skewness, kurtosis, and the scaled-moments-about-the mean up to the K-th. Try:`):
print(`Stat((1+t)^10000,t,8);`):

elif nargs=1 and args[1]=STATg then
print(`STATg(n,a,n1,K): The Statistical information : mean, s.d.,scaled moments up to the K't of the duration of a Markov process given by GFg(m,a,x,n1) (q.v.). Try:`):
print(`STATg(n,n*a^n,10,4);`):

elif nargs=1 and args[1]=STATrn then
print(`STATrn(r,n,K): The EXACT average, s.d., and the third, ..., up to the K2-th moment. Try:`):
print(`STATrn(10,10,6);`):


elif nargs=1 and args[1]=SymMoms then
print(`SymMoms(n,F,a,K): Given a symbol n, and an expression, a (that must be a polynomial in n a^n (a^(n-1)) outputs`):
print(`expilcit expression for the expectation, variance followed by the  moments ABOUT THE MEAN up to the K-th`):
print(`of the r.v. duration of the Markov process of GFg(n,a,x,n1) (q.v.). Try:`):
print(`SymMoms(n,a^n,a,4);`):

elif nargs=1 and args[1]=SymMomsS then
print(`SymMomsS(n,F,a,K): Given a symbol n, and an expression, a (that must be a polynomial in n and a^n (e.g. n*a^(n-1)) outputs`):
print(`expilcit expression for the expectation, variance followed by the SCALED moments (about the mean) up to the K-th`):
pring(`Except that for the odd ones (starting with the third), we give the square, so as not to fuss with the square-root sign.`):
print(`of the r.v. duration of the Markov process of GFg(n,a,x,n1) (q.v.). Try:`):
print(`SymMomsS(n,a^n,a,4);`):


elif nargs=1 and args[1]=SymMomsSlim then
print(`SymMomsSlim(n,F,a,K): Given a symbol n, and an expression F,  in a^n and a  that must be a polynomial in  and a^n (e.g. a^n) outputs`):
print(`expilcit expression for the expectation, variance followed by the LIMITS as n goes to infinity of the scaled moments up to the K-th`):
print(`of the r.v. duration of the Markov process of GFg(n,a,x,n1) (q.v.). Try:`):
print(`SymMomsSlim(n,a^n,a,4);`):

elif nargs=1 and args[1]=TP then
print(`TP(r,n,x): The prob. generating function for going from state (r,n) to state (r-t,n) for t=0..r, using the Dym-Luks expression (that uses inclusion-exclusion).T ry:`):
print(`TP(5,5,x); `):

elif nargs=1 and args[1]=TPbf then
print(`TPbf(r,n,x): The prob. generating function for going from state (r,n) to state (r-t,n) for t=0..r. By Brute Force. For checking only. Try:`):
print(`TPbf(5,5,x); `):


else
 print(`There is no such thing as`, args):

fi:


end:


with(combinat):


#RF(r,n): a random placing of balls {1,..,r} into cells {1,...,n} it outputs a list of length r where L[i] is the cell where ball is placed. Try:
#RF(5,7);
RF:=proc(r,n) local ra,i:
ra:=rand(1..n):
[seq(ra(),i=1..r)]:
end:


#NuCap(L): Given a list of integers L outputs the number of entries that show up exactly once. Try:
#nuCap([5,4,2,2,1,1,3]);
NuCap:=proc(L) local x,co,i,S,gu:

S:={seq(L[i],i=1..nops(L))}:
gu:=add(x[L[i]],i=1..nops(L)):

co:=0:
for i in S do
if coeff(gu,x[i],1)=1 then
 co:=co+1:
fi:
od:
co:

end:

#OSold(r,n): One step of the Dym/Luks game if there are r balls, and n cells. It returns the number of balls not captured
OSold:=proc(r,n) local gu,x,co,i:
gu:=RF(r,n):
gu:=add(x[gu[i]],i=1..r):
co:=0:
for i from 1 to n do
if coeff(gu,x[i],1)=1 then
 co:=co+1:
fi:
od:
r-co:
end:


#OS(r,n): One step of the Dym/Luks game if there are r balls, and n cells. It returns the number of balls not captured
OS:=proc(r,n):r-NuCap(RF(r,n)):end:

#DL(r,n): Playing the Dym-Luks Ball and Cell game and recording its duration. Try:
#DL(10,10);
DL:=proc(r,n) local r1,i:
r1:=r:

for i from 1 while r1>0 do
r1:=OS(r1,n):
od:

i-1:

end:




#OSv(r,n): Verbose version of OS(r,n) (q.v.). Try:
#OSv(100,100);
OSv:=proc(r,n) local L,a,S,gu,i,co:
L:=RF(r,n):

print(`The`, r , `balls landed in the following cells`):
print(``):
print(L):


S:={seq(L[i],i=1..nops(L))}:
gu:=add(x[L[i]],i=1..nops(L)):
print(``):
print(`The set of cells that are occupied is`):
print(``):
print(S):
print(``):
co:={}:
for i in S do
if coeff(gu,x[i],1)=1 then
 co:=co union {i}:
fi:
od:
print(``):
print(`The set of cells that have exactly one ball in them is`):
print(``):
print(co):
print(``):
print(`We remove these balls and are now left with`, r-nops(co), `balls `):
r-nops(co):
end:


#DLv(r,n): A verbose version of DL(r,n) (q.v.). Try:
#DLv(10,10);
DLv:=proc(r,n) local r1,i:
r1:=r:

for i from 1 while r1>0 do
r1:=OSv(r1,n):
od:

print(`This took`, i-1,  `rounds `):

end:


#SimuGrnF(r,n,x,K): an approximation to GFrn(r,n,x) (q.v.) by playing the game K times. Try:
#Of course, every time you get a (slightly) different output (since this is random). As K gets bigger, of course, 
#you should get more agreement. Try:
#SimuGFrn(3,3,x,100);
SimuGFrn:=proc(r,n,x,K) local i:
evalf(add(x^DL(r,n),i=1..K)/K):
end:

#Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try:
#Stat((1+t)^10000,4);
Stat:=proc(F,t,K) local gu,mu,sd,k,f:
if diff(F,t)=0 then
 RETURN([F,0$(K-1)]):
fi:
f:=F/subs(t=1,F):
mu:=subs(t=1,diff(f,t)):
gu:=[mu]:
f:=f/t^mu:
f:=t*diff(t*diff(f,t),t):
sd:=sqrt(subs(t=1,f)):
gu:=[op(gu),sd]:
if sd=0 then
 RETURN(normal([op(gu),0$(K-2)])):
fi:

for k from 3 to K do
 f:=t*diff(f,t):
gu:=expand([op(gu),subs(t=1,f)/sd^k]):
od:

normal(gu):
end:



#SimuSTATrn(r,n,K1,K2): Simulates the Dym-Luks game K1 time and finds the average, s.d., and the third, ..., up to the K2-th moment. Try:
#SimuSTATrn(10,10,100,6);
SimuSTATrn:=proc(r,n,K1,K2) local x:
evalf(Stat(SimuGFrn(r,n,x,K1),x,K2)):
end:

#FUNrn(r,n): all the functions from {1,..,r} to {1,..,n} given as lists of length r with entries in {1,...,n}. Try:
#FUNrn(2,2)
FUNrn:=proc(r,n) local gu,i,gu1:
option remember:
if r=1 then
 RETURN({seq([i],i=1..n)}):
fi:
gu:=FUNrn(r-1,n):
{seq(seq([op(gu1),i],gu1 in gu),i=1..n)}:
end:

#TPbf(r,n,x): The prob. generating function for going from state (r,n) to state (r-t,n) for t=0..r. By Brute Force. For checking only. Try:
#TPbf(5,5,x) 
TPbf:=proc(r,n,x) local gu,gu1:
gu:=FUNrn(r,n):
add(x^(r-NuCap(gu1)),gu1 in gu)/nops(gu):
end:


TP:=proc(r,n,x) local t,j:
expand(add(add((-1)^(j-t)*binomial(j,t)*binomial(n,j)*binomial(r,j)*j!*(n-j)^(r-j)/n^r,j=t..r)*x^(r-t),t=0..r)):
end:


#Prt(n,r,t): The probability of getting from state r to state r-t, what the Dym-Luks paper calls P_{r,r-t}(n). Try: Prt(2,r,1);
Prt:=proc(n,r,t) local j :
expand(add((-1)^(j-t)*binomial(j,t)*binomial(n,j)*binomial(r,j)*j!*(n-j)^(r-j)/n^r,j=t..n)):
end:




#GFrn(r,n,x): the rational function in x giving the probability generating function for the duration of the Dym-Luks Ball and Cell Game with r balls and n cells. Try:
#GFrn(4,4,x);
GFrn:=proc(r,n,x) local gu,t,z,lu:

option remember:

if r=0 then
 RETURN(1):
fi:


gu:=TP(r,n,z):

lu:=0:

for t from 1 to r do lu:=normal(lu+coeff(gu,z,r-t)*GFrn(r-t,n,x)): od:
normal(x*lu/(1-coeff(gu,z,r)*x)):

end:


#GFrnA(r,n,x): the rational function that is the approximation to GFrn(r,n,x) for small n and large r. Easier to deal with. Try:
#GFrnA(4,4,x);
GFrnA:=proc(r,n,x) local lu,gu:
option remember:
if r<=n then
RETURN(GFrn(r,n,x)):
else
gu:=GFrnA(r-1,n,x):
lu:=r*((n-1)/n)^(r-1):
gu:=x*gu*lu/(1-x*(1-lu)):
gu:=normal(gu):
RETURN(gu):
fi:

end:


#GFrn(r,n,x): the rational function in x giving the probability generating function for the duration of the Dym-Luks Ball and Cell Game with r balls and n cells. Try:
#GFrn(4,4,x);
GFrnN:=proc(r,n,x) local gu,t,z:
option remember:

if r=0 then RETURN(1): fi:

gu:=TP(r,n,z):

normal(x*add(coeff(gu,z,r-t)*GFrn(r-t,n,x),t=1..r)/(1-x*coeff(gu,z,r))):


end:



#STATrn(r,n,K): The EXACT average, s.d., and the third, ..., up to the K2-th moment. Try:
#STATrn(10,10,6);
STATrn:=proc(r,n,K) local x:
normal(Stat(GFrn(r,n,x),x,K)):
end:



#STATrnA(r,n,K): The APPROX. average, s.d., and the third, ..., up to the K2-th moment. Using the meomory-one approximation. Try:
#STATrnA(10,10,6);
STATrnA:=proc(r,n,K) local x:
normal(Stat(GFrnA(r,n,x),x,K)):
end:

#PlotDist(r,n,K): plots the distribution of the duration of the Dym-Luks Ball and Cell game with r balls and n cells, using up to the K-th Taylor expansion. Try:
#PlotDist(10,10,30);
PlotDist:=proc(r,n,K) local gu,i,x:
gu:=GFrn(r,n,x):
gu:=taylor(gu,x=0,K+1):
plot([seq([i,evalf(coeff(gu,x,i))],i=1..K)]):
end:


#PlotDist1(r,n,K): plots the distribution of the duration of the Dym-Luks Ball and Cell game with r balls and n cells, using up to the K-th Taylor expansion. Try:
#plotDist1(10,10,30);
PlotDist1:=proc(r,n,K) local gu,i,x:
gu:=GFrn(r,n,x):
gu:=taylor(gu,x=0,K+1):
plot([seq([i,evalf(coeff(gu,x,i))],i=1..K)], scaling=constrained):
end:

#LimDL(n,R,K): For a fixed (numeric) n, and positive integer R, and positive integer K,
# finds K lists. The first list consists of the first R terms M_n(r) (the average duration of the Dym-Luks game with r balls and n cells) for r from 2 to R. 
#The second list of the ratio of s.d/average, and the 3rd through K lists the scaled moments. Try:
#LimDL(3,r,4)
LimDL:=proc(n,R,K) local r,gu,i,lu,j:

for i from 1 to K do
gu[i]:=[]:
od:

for r from 1 to R do
 lu:=STATrn(r,n,K):
 gu[1]:=[op(gu[1]),evalf(lu[1]/((n/r)*(n/(n-1))^(r-1)))]:
 gu[2]:=[op(gu[2]),evalf(lu[2]/lu[1])]:
for j from 3 to K do
 gu[j]:=[op(gu[j]),evalf(lu[j])]:
od:
od:

[seq(gu[j],j=1..K)]:
end:



#Mnr(n,r): STATrn(r,n,2)[1], i.e. the average duration of a Dym-Luks game with  r balls and n cells. Try:
#Mnr(2,10);
Mnr:=proc(n,r): STATrn(r,n,2)[1]:end:


#M2nr(n,r): STATrn(r,n,2)[2]^2, i.e. the average duration of a Dym-Luks game with  r balls and n cells. Try:
#M2nr(2,10);
M2nr:=proc(n,r): STATrn(r,n,2)[2]^2:end:


#Enr(n,r): Mnr(n,r)-Sum((1/j)*(n/(n-1))^(j-1),j=1..r): Try
#Enr(2,10);
Enr:=proc(n,r) local j:Mnr(n,r)-add((n/(n-1))^(j-1)/j,j=1..r):end:


#EnrR(n,r): Mnr(n,r)/Sum((1/j)*(n/(n-1))^(j-1),j=1..r): Try
#EnrR(2,10);
EnrR:=proc(n,r) local j:Mnr(n,r)/add((n/(n-1))^(j-1)/j,j=1..r):end:


#E2nr(n,r): M2nr(n,r)- (Sum(((1/j)*(n/(n-1))^(j-1))^2,j=1..r)-Sum((1/j)*(n/(n-1))^(j-1),j=1..r))
#E2nr(2,10);
E2nr:=proc(n,r) local j:M2nr(n,r)+add((n/(n-1))^(j-1)/j,j=1..r)-add(((n/(n-1))^(j-1)/j)^2,j=1..r):
end:


#E2nrR(n,r): M2nr(n,r)/(Sum(((1/j)*(n/(n-1))^(j-1))^2,j=1..r)-Sum((1/j)*(n/(n-1))^(j-1),j=1..r))
#E2nrR(2,10);
E2nrR:=proc(n,r) local j:M2nr(n,r)/(-add((n/(n-1))^(j-1)/j,j=1..r)+add(((n/(n-1))^(j-1)/j)^2,j=1..r)):
end:


#Paper1(R,n,K): A paper about the Dym-Luks Ball and Cell Game with fixed r and general (symbolic) n. For r=1,...R, it gives the
#prob. generating function and all the moments up to the K's
Paper1:=proc(R,n,K) local lu,gu,x,r,t0:

t0:=time():

print(`The Statistical Distribution of the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed r from 1 to`, R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that you have r balls and (uniformly at) randomly throw them into n cells. Whenever a ball has no cell-mates, it is removed from the game and this continues until all the balls are removed`):
print(``):
print(`We will give the prob. generating function for the duration of the game, i.e. the rational function whose Taylor coefficient of x^i is the probability of the game lasting i rounds is`):
print(``):
print(` for r from 1 to `, R, `as well as the expectation, stanard-deviation, and the scaled moments (about the mean) up to the `, K):
print(``):

for r from 1 to R do

gu:=GFrn(r,n,x):
lu:=STATrn(r,n,K):
print(``):

print(`When there are`, r, `balls then the prob. gen. function, in x,  of the duration is`):
print(``):
print(gu):
print(``):
print(`and in Maple notation:`):
print(``):
lprint(gu):
print(``):
print(`The expected number of rounds is`, lu[1], `the standard-deviation is`, lu[2], `and the scaled moments up to the`, K, `-th are `):
print(``):
print([op(3..K,lu)]):
print(``):
print(`Again in Maple notation,  altogether it is:`):
print(``):
lprint(lu):
print(``):
print(`---------------------------------------------`):
od:

print(`-----------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to produce. `):
print(``):

end:



#Paper2(R,N,K): A paper about the Dym-Luks Ball and Cell Game with fixed (small) n from 1 to N,  and r from 1 to R, it gives the
#prob. generating function and all the moments up to the K's. Try:
#Paper2(5,5,4);
Paper2:=proc(R,N,K) local n,lu,gu,x,r,t0:

t0:=time():

print(`The trends in the Statistical Distribution of the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed n from 1 to`, N, `as r grows larger `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that you have r balls and (uniformly at) randomly throw them into n cells. Whenever a ball has no cell-mates, it is removed from the game and this continues until all the balls are removed`):
print(``):
print(`We will give the prob. generating function for the duration of the game, i.e. the rational function whose Taylor coefficient of x^i is the probability of the game lasting i rounds is`):
print(``):
print(` for n from 2 to`, N, `and for r from 1 to `, R, `as well as the expectation, stanard-deviation, and the scaled moments (about the mean) up to the `, K):
print(``):

for n from 2 to N do
print(`If there are`, n, `cells, then the sequence of rational functions in x giving the prob. generating function of the Dym-Luks game for r  from 1 to`, R, `are as follows `):
print(``):


gu:=[seq(GFrn(r,n,x),r=1..R)]:

print(``):
print(gu):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu):
print(``):
lu:=[seq(evalf(STATrn(r,n,K)),r=1..R)]:
print(`The sequence of the expectation, s.d., followed by the scaled  moments (about the mean) up to the `, K, `are: `):
print(``):
print(lu):
print(``):
print(`The sequence of E_n(r) (the difference with the approximation), is `):
print(``):
lprint([seq(evalf(Enr(n,r)),r=1..R)]):
print(``):
print(`---------------------------------------------`):
od:

print(`-----------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to produce. `):
print(``):

end:





#Paper3(R,a,K): A paper about the Dym-Luks Ball and Cell Game with r balls and trunc(a*r) cells
#prob. generating function and all the moments up to the K's. Try:
#Paper3(10,1,4);
Paper3:=proc(R,a,K) local n,lu,gu,x,r,t0:

t0:=time():

print(`The trends in the Statistical Distribution of the Duration of the Dym-Luks Ball and Cell Game with r balls  and `, a*r , ` cells  as r grows larger `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that you have r balls and (uniformly at) randomly throw them into `, a*r.  `cells. Whenever a ball has no cell-mates, it is removed from the game and this continues until all the balls are removed`):
print(``):
print(`We will give the prob. generating function for the duration of the game, i.e. the rational function whose Taylor coefficient of x^i is the probability of the game lasting i rounds is`):
print(``):

print(` for r from 2 to`, R):
print(``):

for r from 2 to R do
print(`If there are`, trunc(a*r), `cells, and `, r, `balls,then the  rational function in x giving the prob. generating function of the Dym-Luks game is as follows: `):
print(``):


gu:=GFrn(r,trunc(a*r),x):

print(``):
print(gu):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu):
print(``):
print(`The expectation, s.d, etc. are `):
lu:=evalf(STATrn(r,trunc(a*r),K)):
print(``):
print(lu):
print(``):
print(`---------------------------------------------`):
od:

print(`-----------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to produce. `):
print(``):

end:





#Paper1t(R,n,K): An abbreviated version of Paper1(R,n,K) (q.v.) only the average, s.d. and scaled moments.
#Try:
Paper1t:=proc(R,n,K) local lu,gu,x,r,t0:

t0:=time():

if not type(K,integer) and K>=1 then
RETURN(FAIL):
fi:

if K=1 then
print(`The Expectation of  the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed r from 1 to`, R):
elif K=2 then
print(`The Expectation and Variance of  the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed r from 1 to`, R):
else
print(`The Expectation,  Variance, and scaled moments up to the `, K, `-th  of  the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed r from 1 to`, R):
fi:

print(``):
print(`By Shalosh B. Ekhad `):
print(``):
for r from 1 to R do
print(`Suppose that you have`, r , `balls and (uniformly at) randomly throw them into n cells. Whenever a ball has no cell-mates, it is removed from the game and this continues until all the balls are removed`):
print(``):
lu:=STATrn(r,n,K):
print(``):

print(``):

if K>=1 then
print(`The expected number of rounds is`, lu[1]):

print(``):
print(`and in Maple notation`):
print(``):
lprint(lu[1]):
print(``):

fi:

if K>=2 then
print(`the variance is`, normal(lu[2]^2)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu[2]):
print(``):


fi:

if K>=3 then
print(``):
print(`The skewness (aka the scaled third-moment about the mean is`, lu[3]):

print(``):
print(`and in Maple notation`):
print(``):
lprint(lu[3]):
print(``):
fi:


if K>=4 then
print(``):
print(`The kurtosis, aka the scaled fourth-moment about the mean is`, lu[4]):

print(``):
print(`and in Maple notation`):
print(``):
lprint(lu[4]):
print(``):
fi:

if K>=5 then
print(`The remaining scaled moments about the mean from the fifth through the`, K, `- th  are (only in Maple notation)`):
print(``):
lprint([op(5..K,lu)]):
print(``):
fi:
od:

print(`-----------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to produce. `):
print(``):

end:





#Paper2t(R,N,K): Like  Paper2(R,N,K) but only with moments. Try
#Paper2t(5,5,4);
Paper2t:=proc(R,N,K) local n,lu,r,t0,i,j:

t0:=time():

if K=1 then
print(`The trends in the Expectation of the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed n from 1 to`, N, `as r grows larger `):

elif K=2 then
print(`The trends in the Expectation and Standard-Deviation of  the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed n from 1 to`, N, `as r grows larger `):
else
print(`The trends in the Expectation, Variance and moemnts up to the`, K, `of the Duration of the Dym-Luks Ball and Cell Game with r balls  and n cells for fixed n from 1 to`, N, `as r grows larger `):
fi:

print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for n from 2 to N do

print(`Suppose that you have r balls and (uniformly at) randomly throw them into`, n , `cells. Whenever a ball has no cell-mates, it is removed from the game and this continues until all the balls are removed`):
print(``):
print(``):


lu:=[seq(evalf(STATrn(r,n,K)),r=1..R)]:


print(`The sequence of the expectation for balls from r=1 to r= `,R, `are `):
print(``):
print([seq(lu[i][1],i=1..nops(lu))]):
print(``):
print(`The sequence of E_n(r), is `):

print(``):
lprint([seq(evalf(Enr(n,r)),r=1..R)]):

if K>=2 then
print(``):
print(`The sequence of the standard-deviations for number of balls up to `, R, `is ` ):
print(``):
print([seq(lu[i][2],i=1..nops(lu))]):

fi:


for j from 3 to K do
print(``):
print(`The sequence of scaled`, j, ` -th moments are `):
print(``):
print([seq(lu[i][j],i=1..nops(lu))]):



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

od:
print(`-----------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to produce. `):
print(``):

end:




#Paper3t(R,a,K): Like Paper3(R,a,K) (q.v.) but only moments (i.e. no prob. gen. functions). Try:
#Paper3t(10,1,4);
Paper3t:=proc(R,a,K) local lu,r,t0:

t0:=time():

print(`The trends in the Statistical Distribution of the Duration of the Dym-Luks Ball and Cell Game with r balls  and `, a*r , ` cells  as r grows larger `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that you have r balls and (uniformly at) randomly throw them into `, a*r.  `cells. Whenever a ball has no cell-mates, it is removed from the game and this continues until all the balls are removed`):
print(``):
print(`We will give the average, standard deviation and  moments up to the `, K , ` of the duration of the game `):
print(``):

print(` for r from 2 to`, R):
print(``):

for r from 2 to R do

print(`If there are `, r, `balls and `, trunc(a*r), `cells then `):
print(``):
print(`The expectation, s.d, followed by the scaled moments up to the`, K, `-th  are `):
print(``):
lu:=evalf(STATrn(r,trunc(a*r),K)):
print(``):
print(lu):
print(``):
print(`---------------------------------------------`):
od:

print(`-----------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to produce. `):
print(``):

end:


#Paper4(N,R): A table of the quantity Enr(n,r) for n from 2 to N and r from 1 to R. Try:
#Paper4(3,40);
Paper4:=proc(N,R) local n,r,j,t0:
t0:=time():
print(`A table of the difference between the Expected Duration of the Dym-Luks Ball and Cell game with n cells and r balls and`):
print(``):
print(Sum((1/j)*(n/(n-1))^(j-1),j=1..r) ):
print(``):
print(`for n from 2 to`, N, `and r from 1 to`, R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let's call the difference E_n(r)`):
print(``):
for n from 2 to N do

print(`when there are`, n, `cells, the first `, R, `terms of the seqence E_n(r) are`):
print(``):
lprint(evalf([seq(Enr(n,r),r=1..R)])):
print(``):
od:
print(``):
print(`This ends this paper that took`, time()-t0, `seconds to generate . `):
print(``):
end:


#EstC(L,d,r): Given a sequence of numbers L believed to have the asympotic expansion C[0]+C[1]/sqrt(r)+C[2]/r+.. estimates assuming that
#It returns the list [C[0],.., C[d]], looking at the last d+1 terms. It also returns looking at d+1 terms before that.
#it ends in d. Try:
#EstC([seq(5+1/sqrt(i)+3/i,i=1..100)],2,r);
EstC:=proc(L,d,r) local C,i,F,eq,var,r1:

var:={seq(C[i],i=0..d)}:

F:=add(C[i]/sqrt(r)^i,i=0..d):

eq:={seq(L[r1]-subs(r=r1,F),r1=nops(L)-d..nops(L))}:

subs(solve(eq,var),F):

end:


#GFg(n,a,x,n1): The duration probability generating function for the Markov process with prob. of n->n-1 a(n) and n->n (1-a(n)). If currently at n1. Try:
#GFg(n,1/2^n,x,10);
GFg:=proc(n,a,x,n1) local lu:
option remember:

if n1=0 then
 RETURN(1):
else
lu:=subs(n=n1,a):
RETURN(normal(GFg(n,a,x,n1-1)*lu*x/(1-x*(1-lu)))):
fi:
end:

#STATg(n,a,n1,K): The Statistical information : mean, s.d.,scaled moments up to the K't of the duration of a Markov process given by GFg(m,a,x,n1) (q.v.). Try:
#STATg(n,n*a^n,10,4);
STATg:=proc(n,a,n1,K) local x:
Stat(GFg(n,a,x,n1),x,K):
end:


#GuessExpPol1(L,d1,d2,a,n): Given a list of expressions in a and a^n, L, finds a polyomial of degree d1 in a^n and degree d2 in n, let's call it POL such that
#L[n1]=POL(n1,a^n1). Try:
#GuessExpPol1([seq(n1*a^n1,n1=1..10)],1,1,a,n);
GuessExpPol1:=proc(L,d1,d2,a,n) local P,c,eq,var,i,j,n1:

if nops(L)<=(d1+1)*(d2+1)+4 then
 print(`The list`, L, `should be at least of length`, (d1+1)*(d2+1)+4):
 RETURN(FAIL):
fi:

P:=add(add(c[i,j]*(a^n)^i*n^j,j=0..d2),i=0..d1):
var:={seq(seq(c[i,j],j=0..d2),i=0..d1)}:

eq:={seq(subs(n=n1,P)-L[n1],n1=1..(d1+1)*(d2+1)+4)}:

var:=solve(eq,var):

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

P:=subs(var,P):

if {seq(normal(subs(n=n1,P)-L[n1]),n1=(d1+1)*(d2+1)+5..nops(L))}<>{0} then
 RETURN(FAIL):
fi:

P:

end:



#GuessExpPol(L,d1,d2,a,n): Given a list of expressions in a and a^n, L, finds a polyomial of degree<= d1 in a^n and degree<= d2 in n, let's call it POL such that
#L[n1]=POL(n1,a^n1). Try:
#GuessExpPol([seq(n1*a^n1,n1=1..10)],1,1,a,n);
GuessExpPol:=proc(L,d1,d2,a,n) local P,d1a,d2a:

if nops(L)<=(d1+1)*(d2+1)+4 then
 print(`The list`, L, `should be at least of length`, (d1+1)*(d2+1)+4):
 RETURN(FAIL):
fi:

for d1a from 0 to d1 do
 for d2a from 0 to d2 do
 P:=GuessExpPol1(L,d1a,d2a,a,n):
  if P<>FAIL then
    RETURN(P):
  fi:
 od:
od:
FAIL:

end:




#SymMoms(n,F,a,K): Given a symbol n, and an expression F,  in a^n and a  that must be a polynomial in  and a^n (e.g. a^n) outputs
#expilcit expression for the expectation, variance followed by the moments about the mean up to the K-th
#of the r.v. duration of the Markov process of GFg(n,a,x,n1) (q.v.). Try:
#SymMoms(n,a^n,a,4);
SymMoms:=proc(n,F,a,K) local M,L,n1,d1,N,Y,mu,x,k,Y1,gu,i,j:
option remember:

d1:=degree(subs(a^n=N,expand(F)),N):

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

Y:=[]:

M:=[seq(GFg(n,F,x,n1),n1=1..K*d1+6)]:

for k from 1 to K do

M:=[seq(x*diff(M[i],x),i=1..nops(M))]:


L:=normal(subs(x=1,M)):

L:=normal(subs(a=1/a,L)):

mu:=GuessExpPol1(L,k*d1,0,a,n):

mu:=normal(subs(a=1/a,mu)):

Y:=[op(Y),mu]:
od:

mu:=Y[1]:

Y1:=[mu]:

for i from 2 to K do
 gu:=normal(expand((-mu)^i+add(binomial(i,j)*(-mu)^j*Y[i-j],j=0..i-1))):
Y1:=[op(Y1),factor(gu)]:
od:





end:



#SymMomsS(n,F,a,K): Given a symbol n, and an expression F,  in a^n and a  that must be a polynomial in  and a^n (e.g. a^n) outputs
#expilcit expression for the expectation, variance followed by the scaled moments up to the K-th (but the odd ones are squared)
#of the r.v. duration of the Markov process of GFg(n,a,x,n1) (q.v.). Try:
#SymMomsS(n,a^n,a,4);
SymMomsS:=proc(n,F,a,K) local gu,i,lu:
option remember:
gu:=SymMoms(n,F,a,K):
lu:=[gu[1],gu[2]]:

for i from 3 to K do
if i mod 2=1 then
 lu:=[op(lu),factor(normal(gu[i]^2/gu[2]^i))]:
else
 lu:=[op(lu),factor(normal(gu[i]/gu[2]^(i/2)))]:
fi:
od:
lu:
end:






#SymMomsSlim(n,F,a,K): Given a symbol n, and an expression F,  in a^n and a  that must be a polynomial in  and a^n (e.g. a^n) outputs
#expilcit expression for the expectation, variance followed by the LIMITS as n goes to infinity of the scaled moments up to the K-th
#of the r.v. duration of the Markov process of GFg(n,a,x,n1) (q.v.). Try:
#SymMomsSlim(n,a^n,a,4);
SymMomsSlim:=proc(n,F,a,K) local gu,i,lu,ka,ka1,ka2,N:
gu:=SymMomsS(n,F,a,K):

lu:=[gu[1],gu[2]]:


for i from 3 to K do
 ka:=gu[i]:
 ka:=simplify(subs(a=1/a,ka)):

 ka1:=expand(numer(ka)):
 ka1:=subs(a^n=N,ka1):

 ka2:=expand(denom(ka)):
 ka2:=subs(a^n=N,ka2):

 if degree(ka1,N)<>degree(ka2,N) then
  RETURN(FAIL):
  else
  if i mod 2 =1 then
   lu:=[op(lu),factor(sqrt(subs(a=1/a,lcoeff(ka1,N))/subs(a=1/a,lcoeff(ka2,N))))]:
  else
   lu:=[op(lu),factor(subs(a=1/a,lcoeff(ka1,N))/subs(a=1/a,lcoeff(ka2,N)))]:
  fi:
fi:
od:
lu:

end:



#Paper5(K): an article about the expectation, variance, and moments up to the K-th (including the limit as n goes to infinity) for the duration`)
#of a Markov process from n to 0 where the prob. of moving from n to n-1 is a^n, and of staying at n is 1-a^n. Try:
#Paper5(4);
Paper5:=proc(K) local gu1,gu2,gu3,n,a,t0,i:
t0:=time():
gu1:=SymMoms(n,a^n,a,K):
gu2:=SymMomsS(n,a^n,a,K):
gu3:=SymMomsSlim(n,a^n,a,K):
print(``):
print(`The Expected Duration, Variance, and all the moments up to the`,K, `of the Duration of a Markov Process on the Integers from n to 0`):
print(``):
print(`where at each location the particle moves from n to n-1 with probability a^n, and stays where it is with probability 1-a^n`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem : Consider a walker, on the discrete positive line, starting at n, and  goes left, to n-1 with prob. a^n, and remains at n with prob. 1-a^n`):

print(``):
print(`The expectation number of moves until she makes it to 0 is`):
print(``):
print(gu1[1]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu1[1]):
print(``):
print(``):
print(`The variance of the number of moves until she makes it to 0 is`):
print(``):
print(gu1[2]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu1[2]):
print(``):

if K>=3 then

print(`The skewness is`):
print(``):
print(sqrt(gu2[3])):
print(``):


print(`and in Maple notation`):
print(``):
lprint(sqrt(gu2[3])):
print(``):
print(`as n goes to infinity, the limiting skewness is`):
print(``):
print(gu3[3]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu3[3]):
print(``):
print(`Here is a plot from a=0 to a=1`):
print(``):
print(plot(gu3[3],a=0..1));
print(``):

fi:



if K>=4 then

print(`The kurtosis  is`):
print(``):
print(gu2[4]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(gu2[4]):
print(``):

print(`as n goes to infinity, the limiting kurtosis is`):
print(``):
print(gu3[4]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu3[4]):
print(``):

print(`Here is a plot from a=0 to a=1`):
print(``):
print(plot(gu3[4],a=0..1));
print(``):

fi:

for  i from 5 to K do

if i mod 2=1 then



print(`The scaled`, i, `th moment about the mean is is`):
print(``):
print(sqrt(gu2[i])):
print(``):


print(`and in Maple notation`):
print(``):
lprint(sqrt(gu2[i])):
print(``):
print(`as n goes to infinity, the limiting`, i,  `scaled moment about the mean is`):
print(``):
print(gu3[i]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu3[i]):
print(``):
print(`Here is a plot from a=0 to a=1`):
print(``):
print(plot(gu3[i],a=0..1));
print(``):

else


print(`The scaled`, i, `th moment about the mean is is`):
print(``):
print(gu2[i]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(gu2[i]):
print(``):
print(`as n goes to infinity, the limiting`, i, `-th scaled moment about the mean is is`):
print(``):
print(gu3[i]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu3[i]):
print(``):
print(`Here is a plot from a=0 to a=1`):
print(``):
print(plot(gu3[i],a=0..1));
print(``):


fi:


od:

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





#Paper4R(N,R): A table of the quantities EnrR(n,r) and E2nrR(n,R) for n from 2 to N and r from 1 to R. Try:
#Paper4R(3,40);
Paper4R:=proc(N,R) local n,r,j,t0:
t0:=time():
print(`A table of the ratios between the Expected Duration of the Dym-Luks Ball and Cell game with n cells and r balls and`):
print(``):
print(Sum((1/j)*(n/(n-1))^(j-1),j=1..r) ):
print(``):
print(`and the ratios between the variance and `):
print(``):
print(Sum((1/j^2)*(n/(n-1))^(2*j-2),j=1..r) -Sum((1/j)*(n/(n-1))^(j-1),j=1..r) ):
print(``):
print(`for n from 2 to`, N, `and r from 1 to`, R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let's call the ratios for the expectation RatE_n(r)`):
print(``):
print(`Let's call the ratios for the variance RatE2_n(r)`):
for n from 2 to N do

print(`when there are`, n, `cells, the first `, R, ` starting at r=3, terms of the seqence RatE_n(r) are`):
print(``):
lprint(evalf([seq(EnrR(n,r),r=3..R+2)])):
print(``):


print(`Also the terms of terms of the seqence RatE2_n(r) are`):
print(``):
lprint(evalf([seq(E2nrR(n,r),r=3..R+2)])):
print(``):

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