######################################################################
##RandComp.txt: Save this file as  RandComp.txt                      #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read RandComp.txt<Enter>                           #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec. 26, 2016
 
print(`Created:  Dec. 26, 2016`):
print(` This is RandComp.txt  `):
print(`It is one of the packages that accompanies the article `):
print(` Automated Derivation of Limiting Distributions Of Combinatorial Random Variables   Whose   Generating Functions  are Rational`):
print(`by Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.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(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: NC, NuO, ROU `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  Exper, RaC `):
 print(` `):



elif nops([args])=1 and op(1,[args])=AvMomsAFl then
print(`AvMomsAfl(F,t,w,J,n): Like AvMomsA(F,t,w,J,n,beta), but in floating-point.`):
print(`It returns a list of length 2. The first is the asymptotic formula (in floating-point) for the`):
print(`straight enumeration, and the second is a list consisting of the average, the variance, all the way to the J-th moment.`):
print(`Try:`):
print(`AvMomsAfl(1/(1-t-t^2*w),t,w,4,n);`):

elif nops([args])=1 and op(1,[args])=Exper then
print(`Exper(S,n1,N,i,m): inputs a list S, a large positive integer n1, and a large positive integer N, and i, a member of S`):
print(`picks N random compositions of n1 using S, and outputs a list of length m, whose first entry is`):
print(`the empirical average, the second the variance, and the remaining ones the scaled empirical moments.`):
print(`All in floating-point`):
print(`It also returms the exact values, of the average and standard deviation, using the formula`):
print(`Try:`):
print(`Exper([1,2],100,1000,1,4);`):

elif nops([args])=1 and op(1,[args])=NC then
print(`NC(S,n): inputs a set S of positive integer and a non-negative integer, outputs the number of compositions of n`):
print(`into parts using the members of S. Try:`):
print(`NC({1,2},5);`):

elif nops([args])=1 and op(1,[args])=NuO then
print(`NuO(C,i): inputs a composition C and a pos. integer i, outputs the number of occurences of i`):

elif nops([args])=1 and op(1,[args])=RaC then
print(`RaC(S,n,N): a (uniformly at random) random composition of n using the members of S.`):
print(`N is the largest integer needed`):
print(`Try `):
print(`RaC([1,2],1000,1000);`):

elif nops([args])=1 and op(1,[args])=ROU then
print(`ROU(L): inputs a roullete L, picks a location according to the roullete, try:`):
print(`ROU([2,3]);`):


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



###from BiVariateMoms


#FindOpe(L0,L1,n,N, ord1,deg1,MaxSt): inputs two sequences of rational numbers L0,L1, symbols n (the discrete variable), and positive integers,  ord1,
#deg1, and a starting point st1, finds an operator Ope(n,N^(-1)) of degree deg1 in n and order (degree in N^(-1))  ord1 such that
#L1[n]=ope(n,N^(-1))L1[n] for n>=ord1+st1 for st1<=MaxSt. It also returns the starting point.  Try:
#FindOpe([seq(2^i,i=1..10)],[seq(i*2^i-2^(i-1),i=1..10)],n,N,0,1,10);

FindOpe:=proc(L0,L1,n,N, ord1,deg1,MaxSt) local gu,st1:

if nops(L0)<>nops(L1) then
 print(`L0 and L1 should be of the same length`):
 RETURN(FAIL):
fi:

if nops(L0)< (1+ord1)*(1+deg1)+ord1+7+MaxSt then
 print(`The length of L0 (and L1) must be at least `, (1+ord1)*(1+deg1)+ord1+7+MaxSt):
  RETURN(FAIL):
fi:

for st1 from 1 to MaxSt do

 gu:=FindOpe1(L0,L1,n,N, ord1,deg1,st1) :

  if gu<>FAIL then

   RETURN(gu,st1):

 fi:

od:

FAIL:
end:


#FindOpe1(L0,L1,n,N, ord1,deg1,st1): inputs two sequences of rational numbers L0,L1, symbols n (the discrete variable), and positive integers,  ord1,
#deg1, and a starting point st1, finds an operator Ope(n,N^(-1)) of degree deg1 in n and order (degree in N^(-1))  ord1 such that
#L1[n]=ope(n,N^(-1))L1[n] for n>=ord1+st1. Try:
#FindOpe1([seq(2^i,i=1..10)],[seq(i*2^i-2^(i-1),i=1..10)],n,N,0,1,1);

FindOpe1:=proc(L0,L1,n,N, ord1,deg1,st1) local ope,a,eq1,eq,var,n1,i1,j1:

if nops(L0)<>nops(L1) then
 print(`L0 and L1 should be of the same length`):
 RETURN(FAIL):
fi:


if nops(L0)< (1+ord1)*(1+deg1)+ord1+7+st1 then
 print(`The length of L0 (and L1) must be at least `, (1+ord1)*(1+deg1)+ord1+7+st1):
  RETURN(FAIL):
fi:

ope:= add(add(a[i1,j1]*n^i1*N^(-j1),j1=0..ord1),i1=0..deg1):
var:= {seq(seq(a[i1,j1],j1=0..ord1),i1=0..deg1)}:




eq:={}:

for n1 from st1+ord1 to st1+ord1+(1+ord1)*(1+deg1)+6 do

 eq1:=L1[n1]-add(subs(n=n1,coeff(ope,N,-i1))*L0[n1-i1],i1=0..ord1):

 eq:=eq union {eq1}:
od:



var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
ope:=subs(var,ope):
ope:=add(factor(coeff(ope,N,-i1))*N^(-i1),i1=0..ord1):
if {seq(L1[n1]-add(subs(n=n1,coeff(ope,N,-i1))*L0[n1-i1],i1=0..ord1),n1=st1+ord1+(1+ord1)*(1+deg1)+6..nops(L1))}<>{0} then
 RETURN(FAIL):
else
RETURN(ope):
fi:

fi:

end:


#AvMoms(F,t,w,n0,J): inputs a rational function F of variables t, and w, such that the coefficient of t^n is the weight-enumerator, in w,
#according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.
#It also inputs positive integers n0 and J. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]]  (straight moments)
#Try:
#AvMoms(1/(1-t*(1+w)),t,w,10,2);
#AvMoms(MDktw(2,t,w),t,w,10,2);

AvMoms:=proc(F,t,w,n0,J) local F0,L0,L,i,j:
F0:=F:
L0:=[seq(coeff(taylor(subs(w=1,F0),t=0,n0+1),t,i),i=1..n0)]:

L:=[]:

for j from 1 to J do

F0:=w*diff(F0,w):

L:=[op(L), [seq(coeff(taylor(subs(w=1,F0),t=0,n0+1),t,i),i=1..n0)]]:

od:

[L0,L]:

end:

#AsySize(R,t,n,beta): inputs a rational function R of the variable t, that is the  generating function  of some sequence of positive
#integers a(n). Returns an asympototic formula for a(n) in the form A*beta^n, where beta is the largest root of a polynomial in t
#that is also returned (the denominator of R with t->1/t). Try:
#AsySize(1/(1-t-t^2),t,n,beta):
AsySize:=proc(R,t,n,beta) local R1,A,top1,bot1:
R1:=normal(R):
top1:=numer(R1):
bot1:=denom(R1):
A:=-beta*subs(t=1/beta,top1)/subs(t=1/beta,diff(bot1,t)):
A:=normal(A):
bot1:=numer(subs(t=1/t,bot1)):

A:=rem(numer(A),subs(t=beta,bot1),beta)/denom(A):

[A*beta^n,bot1]:
end:

#AvMomsAstraight(F,t,w,J,n,beta): inputs a rational function F of variables t, and w, such that the coefficient of t^n is the weight-enumerator, in w,
#according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.
#It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]] 
#where L0 is the asymptotic formula for the actual size of the n-th set, in terms of the
#growth constant, denoted by  beta, Av, 2ndMom, ..., JthMom are the asymptotic expressions
#in terms of beta  and n for the average, 2ndMom (NOT about the mean), ..., all the way to the J-th moment
#(also not about the mean), followed by the minimal polynomial in t, such that beta is its largest root.
#Try:
#AvMomsAstraight(1/(1-t-t^2*w),t,w,4,n,beta);
AvMomsAstraight:=proc(F,t,w,J,n,beta) local ord1,ku,bot1,mu,lu,n0,N,j,ope1:
ku:=AsySize( normal(subs(w=1,F)),t,n,beta):

bot1:=ku[2]:


ord1:=degree(bot1,t)-1:

n0:=(1+ord1)*(J+1)+ord1+41:


mu:=AvMoms(F,t,w,n0,J):



lu:=[]:

for j from 1 to J do
 ope1:=FindOpe(mu[1],mu[2][j],n,N, ord1,j,30):
  if ope1=FAIL then
    RETURN(FAIL):
  fi:
 ope1:=ope1[1]:

 ope1:=normal(subs(N=beta,ope1)):
 ope1:=rem(numer(ope1),subs(t=beta,bot1),beta)/denom(ope1):

 lu:=[op(lu),ope1]:
od:

[ku[1],lu,bot1,t]:
end:

#AvMomsA(F,t,w,J,n,beta): inputs a rational function F of variables t, and w, such that the coefficient of t^n is the weight-enumerator, in w,
#according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.
#It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]] 
#where L0 is the asymptotic formula for the actual size of the n-th set, in terms of the
#growth constant, denoted by  beta, Av, 2ndMom, ..., JthMom are the asymptotic expressions
#in terms of beta  and n for the average, 2ndMom (about the mean), ..., all the way to the J-th moment
#(also about the mean), followed by the minimal polynomial in t, such that beta is its largest root.
#Try:
#AvMomsA(1/(1-t-t^2*w),t,w,4,n,beta);
AvMomsA:=proc(F,t,w,J,n,beta) local  lu1,lu,mu,ku,r,gu,j,bot1:
option remember:

gu:=AvMomsAstraight(F,t,w,J,n,beta):

bot1:=gu[3]:


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



lu:=gu[2]:

mu:=lu[1]:

lu1:=[mu]:


for j from 2 to J do
 ku:=add(binomial(j,r)*(-mu)^r*lu[j-r],r=0..j-1)+(-mu)^j:
 ku:=expand(ku):
 ku:=normal(subs(N=beta,ku)):
 ku:=normal(rem(numer(ku),subs(t=beta,bot1),beta)/rem(denom(ku), subs(t=beta,bot1),beta)):

 lu1:=[op(lu1),ku]:
od:

[gu[1],lu1,subs(t=beta,gu[3]),beta]:
end:







#ShoreshG(P,t): inputs a polynomial P in t, finds the largest root in floating-point. Try:
#Shoresh(t^3-t-1,t);
ShoreshG:=proc(P,t) local b1,aluf,si,muam,i:
b1:=evalf([solve(P,t)]):

aluf:=b1[1]:
si:=abs(b1[1]):


for i from 2 to nops(b1) do
 muam:=b1[i]:
 if abs(muam)>si then
   aluf:=muam:
   si:=abs(muam):
 fi:
od:

if abs(coeff(aluf,I,1))>10^(-5) then
 print(`largest root does not seem to be real`):
 RETURN(FAIL):
fi:
aluf:=coeff(aluf,I,0):
aluf:

end:

#AvMomsAfl(F,t,w,J,n): Like AvMomsA(F,t,w,J,n,beta), but in floating-point.
#It returns a list of length 2. The first is the asymptotic formula (in floating-point) for the
#straight enumeration, and the second is a list consisting of the average, the variance, all the way to the J-th moment.
#Try:
#AvMomsAfl(1/(1-t-t^2*w),t,w,4,n);
AvMomsAfl:=proc(F,t,w,J,n) local beta,b1,gu,muam,aluf,i,si:
option remember:
gu:=AvMomsA(F,t,w,J,n,beta):

b1:=ShoreshG(gu[3],beta):

[subs(beta=b1,gu[1]),subs(beta=b1,gu[2])]:
end:

###End from BiVariateMoms


#NC(S,n): inputs a set S of positive integer and a non-negative integer, outputs the number of compositions of n
#into parts using the members of S. Try:
#NC({1,2},5);
NC:=proc(S,n)  local s:
option remember:
if n<0 then
0:
elif n=0 then
 1:
else
add(NC(S,n-s),s in S):
fi:
end:





#ROU(L): inputs a roullete L, picks a location according to the roullete, try:
#ROU([2,3]);
ROU:=proc(L) local i,tot,ra,j:
tot:=convert(L,`+`):

ra:=rand(1..tot)():


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





#RaC(S,n): a (uniformly at random) random composition of n using the members of the list  S.
#Try 
#RaC([1,2],1000);
RaC:=proc(S,n) local gu,mu,i,n1:

n1:=n:

mu:=[]:

while n1>0 do
 gu:=[seq(NC(S,n1-S[i]), i=1..nops(S))]:
 i:=ROU(gu):
 mu:=[op(mu),S[i]]:
 n1:=n1-S[i]:
od:
mu:

end:




#NuO(C,i): inputs a composition C and a pos. integer i, outputs the number of occurences of i
NuO:=proc(C,i) local x,i1:
coeff(add(x[C[i1]],i1=1..nops(C)),x[i],1):
end:



#Exper(S,n1,N,i,m): inputs a list S, a large positive integer n, and a large positive integer N, and i, a member of S
#picks N random compositions of n1 using S, and outputs a list of length m, whose first entry is
#the empirical average, the second the empirical standard deviation, and the remaining ones the scaled empirical moments. up to the m-th
#All in floating-point
#Try:
#Exper([1,2],100,1000,1,4);
Exper:=proc(S,n1,N,i,m) local lu,gu,i1,mu,sd1,ku,j,mj,t,w,F,fu,n:

F:=1/(1-add(t^S[i1],i1=1..nops(S))+t^i-w*t^i   ) :

fu:=AvMomsAfl(F,t,w,2,n)[2]:

fu:=subs(n=n1,fu):


lu:=[seq(NC(S,i1),i1=1..n1)]:
gu:=[seq(NuO(RaC(S,n1),i),i1=1..N) ]:

mu:=evalf(convert(gu,`+`)/N):
sd1:= evalf(sqrt(add((gu[i1]-mu)^2,i1=1..N)/N)):


ku:=[mu,sd1]:

for j from 3 to m do
mj:=evalf(add((gu[i1]-mu)^j,i1=1..N)/N):
mj:=evalf(mj/sd1^j):
ku:=[op(ku),mj]:
od:
ku,[fu[1],sqrt(fu[2])]:
end: