######################################################################
##BINnr.txt: Save this file as  BINnr.txt                            #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  BINnr.txt<Enter>                             #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: 
 
print(`Created:  April 2019`):
print(` This is BINnr.txt `):
print(`It is the Maple package that accompanies the article `):
print(` On the Average Maximal Number of Balls in a  Bin Resulting from Throwing r Balls into n Bins T times `):
print(`by Amir Behrouzi-Far and Doron Zeilberger`):

print(`and also available from 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 Findrec procedures type ezraF();, for help with`):
 print(`a specific procedure, type ezraF(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(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: ABc1, enr `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: ABnrMax ,  AsyEstMax, AsyEstMin ,AveM, BinStory, EstMax, EstMin, RecMax, RecMin, SeqAveMax, SeqAveMin  `):
 print(` `):


elif nops([args])=1 and op(1,[args])=ABc1 then
print(`ABc1(n,r): The Amir Behrouzi-Far constant for the constant C1 such that the max number of balls minus T*r/n is asympt. to C1/sqrt(T)`):

elif nops([args])=1 and op(1,[args])=ABc2 then
print(`ABc2(n,r): The Amir Behrouzi-Far constant for the constant C2 such that the max number of balls minus T*r/n is asympt. to C2/sqrt(T)`):

elif nops([args])=1 and op(1,[args])=ABnrMax then
print(`ABnrMax(n,r,T): the full distribution of the r.v. "max number of balls"  minus r*T/n, with throwing r balls into n boxes T times`):

elif nops([args])=1 and op(1,[args])=ABnrMin then
print(`ABnrMin(n,r,T): the full distribution of the r.v. "min number of balls"  minus r*T/n, with throwing r balls into n boxes T times`):




elif nops([args])=1 and op(1,[args])=AsyEstMax then
print(`AsyEstMax(n,r,MaxC,K1,k): Numerical estimate for the constant C1 such that the expected number of the MAXIMUM number of balls`):
print(`in a bin upon throwing r balls into n bins T times, minus T*r/n is C1*sqrt(T). It uses K1 terms of the sequence`):
print(`obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.`):
print(`It returns the estimate for  K1 in floating-point using an asymptotic solution to order k. Try`):
print(`AsyEstMax(3,1,10,5000,3);`):

elif nops([args])=1 and op(1,[args])=AsyEstMin then
print(`AsyEstMin(n,r,MaxC,K1,k): Numerical estimate for the constant C1 such that the expected number of the MINIMUM number of balls`):
print(`in a bin upon throwing r balls into n bins T times minus T*r/n is C1*sqrt(T). It uses K1 terms of the sequence`):
print(`obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.`):
print(`It returns the estimate for K1 in floating-point using an asymptotic solution to order k. Try`):
print(`AsyEstMin(3,1,10,5000,3);`):

elif nops([args])=1 and op(1,[args])=EstMax then
print(`EstMax(n,r,MaxC,K1): Numerical estimate for the constant C1 such that the expected number of the MAXIMUM number of balls`):
print(`in a bin upon throwing r balls into n bins T times minus T*r/n is C1*sqrt(T). It uses K terms of the sequence`):
print(`of averages to guess a recurrence, and then the K2-th term of the sequence. It uses the K1-th term of the sequence`):
print(`obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.`):
print(`It returns the estimate for  K1 in floating point, just using the sequence. Try`):
print(`EstMax(3,1,10,5000);`):

elif nops([args])=1 and op(1,[args])=EstMin then
print(`EstMin(n,r,MaxC,K1): Numerical estimate for the constant C1 such that the expected number of the MINIMUM number of balls`):
print(`in a bin upon throwing r balls into n bins T times minus T*r/n is C1*sqrt(T). It uses K terms of the sequence`):
print(`of averages to guess a recurrence, and then the K2-th term of the sequence. It uses the K1-th term of the sequence`):
print(`obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.`):
print(`It returns the estimate for K1-1000 and K1 in floating point, just using the sequence. Try`):
print(`EstMin(3,1,10,5000);`):

elif nops([args])=1 and op(1,[args])=AveM then
print(`AveM(L,m): inputs a discrete prob. distribution, outputs the average, variance, and scaled moments`):

elif nops([args])=1 and op(1,[args])=BinStory then
print(`BinStory(n,r,MaxC,K1): tells a story about the average of the largest number of balls in a bin upon throwing`):
print(`r balls into n bins T times. Try:`):
print(`BinStory(3,1,10,5000);`):

elif nops([args])=1 and op(1,[args])=enr then
print(`enr(n,r,x): the elementary symmetric function e_r(x[1], ..., x[n]). Try: `):
print(`enr(4,2,x); `):
print(``):

elif nops([args])=1 and op(1,[args])=RecMax then
print(`RecMax(n,r,T,S,MaxC): guesses a recurrence operator of complexity <=MaxC for the sequence of averages for the r.v. "MAXIMUM number of balls" upon throwing`):
print(`r balls (w/o replacements) into n bins, T times. It uses K data points, and S is the shift operator in T.`):
print(`It returns FAIL if K is not big enough. It also returns the list of initial values needed. Try:`):
print(`RecMax(3,1,T,S,10);`):


elif nops([args])=1 and op(1,[args])=RecMin then
print(`RecMin(n,r,T,S,MaxC): guesses a recurrence operator of complexity <=MaxC for the sequence of averages for the r.v. "MINIMUM number of balls" upon throwing`):
print(`r balls (w/o replacements) into n bins, T times. It uses K data points, and S is the shift operator in T.`):
print(`It returns FAIL if K is not big enough.  It also returns the list of initial values needed. Try:`):
print(`RecMin(3,1,T,S,10);`):


elif nops([args])=1 and op(1,[args])=SeqAveMax then
print(`SeqAveMax(n,r,K): The sequence of averages for the maximum number of balls in a bin when your throw r balls`):
 print(`into n bins (without replacement), T times , minus r*T/n, from T=1 toT= K. Try:`):
print(`SeqAveMax(3,1,100);`):


elif nops([args])=1 and op(1,[args])=SeqAveMin then
print(`SeqAveMin(n,r,K): The sequence of averages for the minimum number of balls in a bin when your throw r balls`):
 print(`into n bins (without replacement), T times minus r*T/n, , from T=1 toT= K. Try:`):
print(`SeqAveMax(3,1,100);`):


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



ez:=proc(): print(` ABc1(n,r), ABc2(n,r)`): end:

#from Findrec 
ezraF:=proc()
if args=NULL then

print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF, SeqFromRec, SeqFromRecFw `):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,Ini,n,N,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,[1],n,N,10);`):


elif nargs=1 and args[1]=SeqFromRecFw then
print(`SeqFromRecFw(ope,Ini,n,N,K,w): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends outputs the values from K-w+1 to K`):
print(`For example, try:`):
print(`SeqFromRecFw(N-n-1,[1],n,N,10,2);`):

fi:
 
end:


###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:
option remember:
for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,Ini,n,N,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,Ini,n,N,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:

 
#SeqFromRecFw(ope,Ini,n,N,K,w): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#outputs the values from K to K+w-1
SeqFromRecFw:=proc(ope,Ini,n,N,K,w)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=evalf(Ini):

for n1 from nops(Ini)+1 to w do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),j1=1..L)]:
od:

for n1 from w+1 to K+w-1 do
 gu:=evalf([op(2..nops(gu),gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),j1=1..L)]):
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
#end from Findrec 

enr:=proc(n,r,x) local t,i:
expand(coeff(mul(1+x[i]*t,i=1..n),t,r)):
end:

#ABnrMax(n,r,T): the full distribution of the r.v. "max number of balls" minus r*T/n, with throwing r balls into n boxes T times
ABnrMax:=proc(n,r,T) local gu,x,T1,gu1,mu,i,i1:
gu:=expand(enr(n,r,x)^T/binomial(n,r)^T ):

for i from trunc(T*r/n) to T do
T1[i]:=0:
od:

for i from 1 to nops(gu) do
gu1:=op(i,gu):
mu:=max(seq(degree(gu1,x[i1]),i1=1..n)):
T1[mu]:=T1[mu]+op(1,gu1):
od:

[seq( [i-T*r/n,T1[i]] ,i=trunc(T*r/n)..T )]:

end:


#ABnrMin(n,r,T): the full distribution of the r.v. "max number of balls" with throwing r balls into n boxes T times
ABnrMin:=proc(n,r,T) local gu,x,T1,gu1,mu,i,i1:
gu:=expand(enr(n,r,x)^T/binomial(n,r)^T ):

for i from 0 to T do
T1[i]:=0:
od:

for i from 1 to nops(gu) do
gu1:=op(i,gu):
mu:=min(seq(degree(gu1,x[i1]),i1=1..n)):
T1[mu]:=T1[mu]+op(1,gu1):
od:

[seq( [i-r*T/n,T1[i]] ,i=0..T )]:

end:

#AveM(L,m): inputs a discrete prob. distribution, outputs the average, variance, and scaled moments
AveM:=proc(L,m) local mu,sig2,i,gu,i1:
if add(L[i][2],i=1..nops(L))<>1 then
 RETURN(FAIL):
fi:


mu:=add(L[i][1]*L[i][2],i=1..nops(L)):

if m=1 then
 RETURN([mu]):
fi:

sig2:=add((L[i][1]-mu)^2*L[i][2],i=1..nops(L)):

gu:=[mu,sig2]:

for i from 3 to m do
 gu:=[op(gu),evalf(add((L[i1][1]-mu)^i*L[i][2],i1=1..nops(L)))/evalf(sig2^(i/2))]:
od:

gu:

end:

#SeqAveMax(n,r,K): The sequence of averages for the maximum number of balls from 1 to K. Try:
#SeqAveNax(3,1,100);
SeqAveMax:=proc(n,r,K) local i:
[seq(AveM(ABnrMax(n,r,i),1)[1],i=1..K)]:
end:

#SeqAveMin(n,r,K): The sequence of averages for the maximum number of balls from 1 to K
SeqAveMin:=proc(n,r,K) local i:
[seq(AveM(ABnrMin(n,r,i),1)[1],i=1..K)]:
end:


 
#ABc1(n,r): The Amir Behrouzi-Far constant for the max. number of balls with r balls and n bins
ABc1:=proc(n,r)
(r/n)*sqrt(Pi*log(n))*log(n/r):
#r/n*sqrt((3*r)^(1/2)*log(n)/n^(1/4))*log(n/r):
end:





#RecMax(n,r,T,S,MaxC): guesses a recurrence operator of complexity <=MaxC for the sequence of averages for the r.v. "max number of balls" upon throwing
#r balls (w/o replacements) into n bins, T times. It uses K data points, and S is the shift operator in T.
#It returns FAIL if K is not big enough. Try:
#RecMax(3,1,T,S,12);
RecMax:=proc(n,r,T,S,MaxC) local i,gu,K,ope:
option remember:
K:=max(seq((2+i)*(1+MaxC-i)+5,i=0..MaxC)):
gu:=SeqAveMax(n,r,K):
ope:=Findrec(gu,T,S,MaxC):

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

ope,SeqAveMax(n,r,degree(ope,S)):

end:


#RecMin(n,r,T,S,MaxC): guesses a recurrence operator of complexity <=MaxC for the sequence of averages for the r.v. "min number of balls" upon throwing
#r balls (w/o replacements) into n bins, T times. It uses K data points, and S is the shift operator in T.
#It returns FAIL if K is not big enough. Try:
#RecMin(3,1,T,S,12);
RecMin:=proc(n,r,T,S,MaxC) local i,gu,K,ope:
option remember:
K:=max(seq((2+i)*(1+MaxC-i)+5,i=0..MaxC)):
gu:=SeqAveMin(n,r,K):
ope:=Findrec(gu,T,S,MaxC):


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

ope,SeqAveMax(n,r,degree(ope,S)):

end:


#ForSol(ope,n,N,k): a formal solution to the recurrence starting  with n^(1/2) to k terms
ForSol:=proc(ope,n,N,k) local c,f,i,gu,	 x,eq,var:

f:=1+add(c[i]/n^i,i=1..k):

gu:=add(normal(subs(n=1/x,coeff(ope,N,i))*subs(n=1/x,subs(n=n+i,f)) )*sqrt(1+i*x),i=0..degree(ope,N)):
gu:=taylor(gu,x=0,2*k+1):
eq:={seq(coeff(gu,x,i),i=k+1..2*k)}:
var:={seq(c[i],i=1..k)}:
var:=solve(eq,var);
subs(var,f):
end:

#AsyEstMax(n,r,MaxC,K1,k): Numerical estimate for the constant C1 such that the expected number of the maximum number of balls
#in a bin upon throwing r balls into n bins T times is T*r/n*(1+C1/sqrt(T)). It uses K terms of the sequence
#of averages to guess a recurrence, and then the K2-th term of the sequence. It uses the K1-th term of the sequence
#obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.
#It returns the estimate  K1 in floating-point using an asymptotic solution to order k. Try
#AsyEstMax(3,1,12,5000,3);
AsyEstMax:=proc(n,r,MaxC,K1,k) local ope,T,S,INI,gu,f,lu,K2:

if K1<=2000 then
 print(`K1>2000`):
 RETURN(FAIL):
fi:

ope:=RecMax(n,r,T,S,MaxC)[1]:

f:=ForSol(ope,T,S,k):

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

INI:=SeqAveMax(n,r,degree(ope,S)):
gu:=SeqFromRec(ope,INI,T,S,K1)[-1]:

evalf(gu/sqrt(K1)/subs(T=K1,f)):

end:

#AsyEstMin(n,r,MaxC,K1,k): Numerical estimate for the constant C1 such that the expected number of the minimum number of balls
#in a bin upon throwing r balls into n bins T times is T*r/n*(1+C1/sqrt(T)). It uses K terms of the sequence
#of averages to guess a recurrence, and then the K2-th term of the sequence. It uses the K1-th term of the sequence
#obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.
#It returns the estimate for K1-1000 and K1 in floating-point using an asymptotic solution to order k. Try
#AsyEstMin(3,1,12,5000,3);
AsyEstMin:=proc(n,r,MaxC,K1,k) local ope,T,S,INI,gu,f,lu,K2:

if K1<=2000 then
 print(`K1>2000`):
 RETURN(FAIL):
fi:

ope:=RecMin(n,r,T,S,MaxC)[1]:

f:=ForSol(ope,T,S,k):

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


INI:=SeqAveMin(n,r,degree(ope,S)):
gu:=SeqFromRec(ope,INI,T,S,K1)[-1]:

evalf(gu/sqrt(K1)/subs(T=K1,f)):

end:




#EstMax(n,r,MaxC,K1): Numerical estimate for the constant C1 such that the expected number of the maximum number of balls
#in a bin upon throwing r balls into n bins T times  minus n*r/T is C1/sqrt(T). It uses K terms of the sequence
#of averages to guess a recurrence, and then the K2-th term of the sequence. It uses the K1-th term of the sequence
#obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.
#It returns the estimate for K1-1000 and K1 in floating-point just using the sequence. Try
#EstMax(3,1,12,5000,3);
EstMax:=proc(n,r,MaxC,K1) local ope,T,S,INI,gu,lu,K2:

if K1<=2000 then
 print(`K1>2000`):
 RETURN(FAIL):
fi:

ope:=RecMax(n,r,T,S,MaxC)[1]:


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


INI:=SeqAveMax(n,r,degree(ope,S)):
gu:=SeqFromRec(ope,INI,T,S,K1)[-1]:

evalf(gu/sqrt(K1)):

end:

#EstMin(n,r,MaxC,K1): Numerical estimate for the constant C1 such that the expected number of the minimun number of balls
#in a bin upon throwing r balls into n bins T times  minus n*r/T is C1/sqrt(T). It uses K terms of the sequence
#of averages to guess a recurrence, and then the K2-th term of the sequence. It uses the K1-th term of the sequence
#obtained from the guessed recurrence. If no recurrence of complexity <=MaxC is found it returns FAIL.
#It returns the estimate for K1-1000 and K1 in floating-point just using the sequence. Try
#EstMin(3,1,12,5000,3);
EstMin:=proc(n,r,MaxC,K1) local ope,T,S,INI,gu,lu,K2:

if K1<=2000 then
 print(`K1>2000`):
 RETURN(FAIL):
fi:

ope:=RecMin(n,r,T,S,MaxC)[1]:


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


INI:=SeqAveMin(n,r,degree(ope,S)):


INI:=SeqAveMin(n,r,degree(ope,S)):
gu:=SeqFromRec(ope,INI,T,S,K1)[-1]:

evalf(gu/sqrt(K1)):

end:

#BinStory(n,r,MaxC,K1): tells a story about the average of the largest number of balls in a bin upon throwing
#r balls into n bins T times. Try:
#BinStory(3,1,12,5000);

BinStory:=proc(n,r,MaxC,K1) local ope,T,S,A,ORD,D1,C,i,gu1,gu2,INI:

ope:=RecMax(n,r,T,S,MaxC)[1]:

INI:=SeqAveMax(n,r,degree(ope,S)):

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

ORD:=degree(ope,S):

ope:=expand(subs(T=T-ORD,ope)/S^ORD):
ope:=add(factor(coeff(ope,S,-i))*S^(-i),i=0..ORD):
if coeff(ope,S,0)<>1 then
 RETURN(FAIL):
fi:
ope:=1-ope:


print(`On the average of the number of the Maximal (and Minimal) Number of Balls upon throwing`, r , `balls into `, n, `bins T times `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose you throw, uniformly at random,`, r, `balls into `, n, `boxes  , T  times `):
print(`Let `, A(T), `be the average of the Maximum number of balls minus `, r*T/n):
print(``):
print(`Then A(T) satisfies the following linear recurrence equation with polynomial coefficients`):
print(``):
print(A(T)=add(coeff(ope,S,-i)*A(T-i),i=1..ORD)):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(A(i)=INI[i],i=1..nops(INI))):

print(``):
print(`and in Maple notation`):
print(``):
lprint(A(T)=add(coeff(ope,S,-i)*A(T-i),i=1..ORD)):
print(``):

print(`Using this recurrence we can compute many terms.`):
print(``):

gu1:=EstMax(n,r,MaxC,K1):

gu2:=EstMax(n,r,MaxC,K1+1000):

D1:=trunc(evalf(-log[10](abs(gu1-gu2))))-1:

if type(D1,posint) then
print(`This enables us to estimate that A(T) is asympotically`, C*sqrt(T) ):
print(``):


print(`where C is approximately`, evalf(gu2,D1) ):
fi:

ope:=RecMin(n,r,T,S,MaxC)[1]:
INI:=SeqAveMin(n,r,degree(ope,S)):

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

ORD:=degree(ope,S):

ope:=expand(subs(T=T-ORD,ope)/S^ORD):
ope:=add(factor(coeff(ope,S,-i))*S^(-i),i=0..ORD):
if coeff(ope,S,0)<>1 then
 RETURN(FAIL):
fi:
ope:=1-ope:

print(``):
print(`Suppose you throw, uniformly at random,`, r, `balls into `, n, `boxes  , T  times `):
print(`Let `, A(T), `be the average of the Minimum number of balls minus `, r*T/n):
print(``):
print(`Then A(T) satisfies the following linear recurrence equation with polynomial coefficients`):
print(``):
print(A(T)=add(coeff(ope,S,-i)*A(T-i),i=1..ORD)):
print(``):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(A(i)=INI[i],i=1..nops(INI))):

print(`and in Maple notation`):
print(``):
lprint(A(T)=add(coeff(ope,S,-i)*A(T-i),i=1..ORD)):
print(``):

print(`Using this recurrence we can compute many terms.`):
print(``):


gu1:=EstMin(n,r,MaxC,K1):

gu2:=EstMin(n,r,MaxC,K1+1000):

D1:=trunc(-log[10](abs(gu1-gu2)))-1:

if type(D1,posint) then

print(`This enables us to estimate that A(T) is asympotically`, C*sqrt(T) ):
print(``):
print(`where C is approximately`, evalf(gu2,D1) ):
fi:


end: