######################################################################
##CLT: Save this file as  CLT                                       #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read CLT<Enter>                                    #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Nov. 4, 2008
 
print(`Created:  Nov. 4, 2008`):
print(` This is CLT `):
print(`to (automatically!) discover (and prove!) the asymptotic expressions`):
print(`for the general moments of n repetitions of`):
print(`any prob. distibution given in terms of its moments`):
print(`Using the the polynomial ansatz`):
print(`It is one of the two packages that accompany the article `):
print(`The Automatic Central Limit Theorems Generator (and Much More!) `):
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(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
print(`For help with procedures dealing with specific discrete`):
print(`prob. dist., type ezraS();`):
print(`For supporting procedures, type, ezra1(); with`):
print(`help on any of the procedures listed by ezraS(); and ezra1();`):
print(`type: ezraS(procedure_name);`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(`AsyNorMomEviaFa, AsyNorMomOviaFa,`):
 print( ` Compare, GuessFactor, GuessPol, ProveAsyFaMomE, `):
 print(` ProveAsyFaMomO, Stir2 `):

else
ezra(args):
fi:

end:


ezra:=proc()

if args=NULL then
 print(`The procedures about general prob. dist. are: `):
 print(`AsyFaMomEg, AsyFaMomOg, AsyMomE1, AsyMomEg, `):
 print(`AsyMomEgM,AsyMomOg,AsyMomOgM, AsyMomO1`):
 print(`AsyNorMomE1,AsyNorMomO1,`):
 print(` AsyNorFaMomEg, AsyNorFaMomOg,AsyNorMomEg, AsyNorMomOg`):
 print(`AsyNorMomEgM, AsyNorMomOgM`):
 print(` FactorialInTermsOfM, FactorialMomsG`):
 print(`NorFaMomsEg,NorFaMomsOg, `):
 print(`  PgfToF , ProveAsyFaMomEg, ProveAsyFaMomOg, ProveAsyFaMomG`):
 print(` SipurG .`):


elif nops([args])=1 and op(1,[args])=AsyFaMomEg then
print(`AsyFaMomEg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  Factorial (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbial die whose factorial moments are P[2],P[3], ...`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyFaMomEg(P,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyFaMomOg then
print(`AsyFaMomOg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  Factorial (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbial die whose factorial moments are P[2],P[3], ...`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyFaMomOg(P,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomEg then
print(`AsyMomEg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbial die whose  p.d. has factorial moments P[2],P[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomEg(P,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomEgM then
print(`AsyMomEgM(M,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbial die whose  p.d. has moments M[1]=0,M[2],M[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomEgM(M,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomE1 then
print(`AsyMomE1(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomE1(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomO1 then
print(`AsyMomO1(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomO1(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomOg then
print(`AsyMomOg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbial die whose  p.d. has factorial moments P[2],P[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomOg(P,n,r,3);`):


elif nops([args])=1 and op(1,[args])=AsyMomOgM then
print(`AsyMomOgM(M,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbial die whose  p.d. has moments M[1]=0,M[2],M[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomOgM(M,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorFaMomEg then
print(`AsyNorFaMomEg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r)-th Factorial moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbail die whose mean in 0 and whose`):
print(`factorial moments are P[2],P[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorFaMomEg(P,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorFaMomOg then
print(`AsyNorFaMomOg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r+1)-th Factorial moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a proverbail die whose mean in 0 and whose`):
print(`factorial moments are P[2],P[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorFaMomOg(P,n,r,3);`):


elif nops([args])=1 and op(1,[args])=AsyNorMomE1 then
print(`AsyNorMomE1(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomE1(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomEg then
print(`AsyNorMomEg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r)-th moment`):
print(`of the r.v. obtained by repeating n times `):
print(` a p.d. whose factorial moments are P[2],P[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomEg(P,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomEgM then
print(`AsyNorMomEgM(M,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r)-th moment`):
print(`of the r.v. obtained by repeating n times `):
print(` a p.d. whose  moments are M[1]=0, M[2],M[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomEgM(M,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomO1 then
print(`AsyNorMomO1(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomO1(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomOg then
print(`AsyNorMomOg(P,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times `):
print(` a p.d. whose factorial moments are P[2],P[3], ...`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomOg(P,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomOgM then
print(`AsyNorMomOgM(M,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times `):
print(` a p.d. whose  moments are M[1]=0, M[2],M[3], ...`):
print(`Here n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomOgM(M,n,r,3);`):

elif nops([args])=1 and op(1,[args])=FactorialInTermsOfM then
print(`FactorialInTermsOfM(M,L): Expressing the first L factorial`):
print(`moments in terms of the moments M[1],M[2], ...`):
print(`For example, try: `):
print(`FactorialInTermsOfM(M,5);`):

elif nops([args])=1 and op(1,[args])=FactorialMomsG then
print(`FactorialMomsG(P,n,L): the first L Factorial`):
print(`moments (about the mean)`):
print(`of the distribution obtained by n throws of a proverbail die`):
print(`whose factorial moments are P[2],P[3], ...`):
print(`For example, for the first 10 moments `):
print(`FactorialMomsG(P,n,10);`):


elif nops([args])=1 and op(1,[args])=NorFaMomsEg then
print(`NorFaMomsEg(P,n,L): the first L normalized `):
print(`even Factorial moments (about the mean)`):
print(`i.e. the list whose r-th entry is the (2r)-th Factorial moment`):
print(`divided by the r-th power of the second moment`):
print(`divided by (2r)!/2^r/r!`):
print(`of the distribution obtained by n throws of a `):
print(`proverbial die`):
print(`whose factorial moments are P[2],P[3] (given symbolically)`):
print(`For example, try `):
print(`NorFaMomsEg(P,n,10)`):

elif nops([args])=1 and op(1,[args])=NorFaMomsOg then
print(`NorFaMomsOg(P,n,L): the first L normalized `):
print(`even Factorial moments (about the mean)`):
print(`i.e. the list whose r-th entry is the (2r)-th Factorial moment`):
print(`divided by the r-th power of the second moment`):
print(`divided by (2r)!/2^r/r!`):
print(`of the distribution obtained by n throws of a `):
print(`proverbial die`):
print(`whose factorial moments are P[2],P[3] (given symbolically)`):
print(`For example, try `):
print(`NorFaMomsOg(P,n,10)`):


elif nops([args])=1 and op(1,[args])=PgfToF then
print(`PgfToF(f,x,L): the first L factorial moments`):
print(`(about the mean)`):
print(`of a discrete p.d. whose pgf is f(x).`):
print(`For example, try:`):
print(`PgfToF(p*x+1-p,x,10);`):

elif nops([args])=1 and op(1,[args])=ProveAsyFaMomEg then
print(`ProveAsyFaMomEg(s): Rigorously proves the`):
print(`conjectured asymptotic expression given by AsyFaMomEg(P,n,r,s),`):
print(`i.e., the asymptotic expression of order s for the even`):
print(`factorial moments of the random variable "total score" `):
print(`upon repeating n times a prob. dist. whose factorial`):
print(`moments are P[1]=0, P[2], ... For example, try:`):
print(`ProveAsyFaMomEg(2);`):

elif nops([args])=1 and op(1,[args])=ProveAsyFaMomOg then
print(`ProveAsyFaMomOg(s): Rigorously proves the`):
print(`conjectured asymptotic expression given by AsyFaMomOg(P,n,r,s),`):
print(`i.e., the asymptotic expression of order s for the odd`):
print(`factorial moments of the random variable "total score" `):
print(`upon repeating n times a prob. dist. whose factorial`):
print(`moments are P[1]=0, P[2], ... For example, try:`):
print(`ProveAsyFaMomOg(2);`):

elif nops([args])=1 and op(1,[args])=ProveAsyFaMomG then
print(`ProveAsyFaMomG(s): Rigorously proves the`):
print(`conjectured asymptotic expression given by AsyFaMomOg(P,n,r,s),`):
print(`and AsyFaMomOg(P,n,r,s),`):
print(`i.e., the asymptotic expression of order s for all the`):
print(`factorial moments of the random variable "total score" `):
print(`upon repeating n times a prob. dist. whose factorial`):
print(`moments are P[1]=0, P[2], ... For example, try:`):
print(`ProveAsyFaMomG(2);`):

elif nops([args])=1 and op(1,[args])=SipurG then
print(`SipurG(M,n,r,s): Proves a refined version of the`):
print(`Central Limit Theorem, complete with asymptotics`);
print(`and the fact that if the first few moments`):
print(`coincide with the momenets of the Normal Distribution`):
print(`then the convergence to the Normal Distribution is`):
print(`faster; Do:`):
print(`SipurG(M,n,r,3);`):

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




ezraS:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` AppxProb, Prob, AsyFaMomE, AsyFaMomO`):
 print(` AsyMomE, AsyMomO,AsyNorFaMomE, AsyNorFaMomO`):
 print(`AsyNorMomE, AsyNorMomO, `):
 print(` CorrectedGaussian, FactorialMoms, Moms, NorFaMomsE, NorFaMomsO `):
 print(`NorMomsE, NorMomsO, ProveAsyFaMom`):
 print(`ProveAsyFaMomEverbose,ProveAsyFaMomOverbose`):



elif nops([args])=1 and op(1,[args])=AppxProb then
print(`AppxProb(f,x,n,a,b,s): Given a die whose p.g.f. is f(x)`):
print(`and an integer n, the approximate probability that the sum`):
print(`of the throws is between a and b inclusive.`):
print(`using a corrected Gaussian of order s.`):
print(`(s=0 corresponds to the usual  Normal Appx.)`):
print(`For example, for the appx. prob. of throwing a `):
print(`fair coin 200 times and getting between 90 and 110`):
print(`Heads, using a corrected Gaussian of order 1, type:`):
print(` AppxProb((1+x)/2,x,200,90,110,1);`):


elif nops([args])=1 and op(1,[args])=AsyFaMomE then
print(`AsyFaMomE(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  Factorial (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyFaMomE(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyFaMomO then
print(`AsyFaMomO(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r+1)-th Factorial moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyFaMomO(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomE then
print(`AsyMomE(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomE(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyMomO then
print(`AsyMomO(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the  (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyMomO(1-p+p*x,x,n,r,3);`):



elif nops([args])=1 and op(1,[args])=AsyNorFaMomE then
print(`AsyNorFaMomE(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r)-th Factorial moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorFaMomE(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorFaMomO then
print(`AsyNorFaMomO(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r+1)-th Factorial moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorFaMomO(1-p+p*x,x,n,r,3);`):


elif nops([args])=1 and op(1,[args])=AsyNorMomE then
print(`AsyNorMomE(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomE(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomEviaFa then
print(`AsyNorMomEviaFa(f,x,n,r,s): like AsyNorMomE(f,x,n,r,s)`):
print(`but via the Factorial moments, for which we have a rigorous proof`):
print(`For example, try: AsyNorMomEviaFa(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomO then
print(`AsyNorMomO(f,x,n,r,s): The first s+1 terms in the`):
print(`asymptotic expansion of the normalized (2r+1)-th moment`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,n,r are symbolic, but s is a numeric non-neg. integer. `):
print(`For example, try: AsyNorMomO(1-p+p*x,x,n,r,3);`):

elif nops([args])=1 and op(1,[args])=AsyNorMomOviaFa then
print(`AsyNorMomOviaFa(f,x,n,r,s): like AsyNorMomO(f,x,n,r,s)`):
print(`but via the Factorial moments, for which we have a rigorous proof`):
print(`For example, try: AsyNorMomOviaFa(1-p+p*x,x,n,r,3);`):


elif nops([args])=1 and op(1,[args])=Compare then
print(`Compare(f,x,n1,a,b,s): Compares the errors in the approximations`):
print(`by the corrected Gaussians to the uncorrected ones`):
print(`for exaple try:`):
print(`Compare((1+x)/2,x,100,40,60,3);`):

elif nops([args])=1 and op(1,[args])=CorrectedGaussian then
print(`CorrectedGaussian(f,x,n,s): The `):
print(`Corrected Gaussian to the limiting distribution`):
print(`of a symmetric p.d. using the`):
print(`asymptotic expansion of the even moments`):
print(`the odd moments are 0 since the dist. is symmetric`):
print(`of the r.v. obtained by repeating n times throwing`):
print(` a die whose probability generating function is f(x)`):
print(`Here x,and n are symbolic, but s is a numeric non-neg. integer `):
print(`For example, try: CorrectedGaussian((1+x)/2,x,n,3);`):


elif nops([args])=1 and op(1,[args])=FactorialMoms then
print(`FactorialMoms(f,x,n,L): the first L Factorial`):
print(`moments (about the mean)`):
print(`of the distribution obtained by n throws of a die`):
print(`whose prob. generating function is f(x)`):
print(`For example, for the first 10 moments `):
print(`for tossing a loaded coin whose`):
print(` probability of Head is p, and the r.v. is [number of Heads], type: ` ):
print(`FactorialMoms(p*x+1-p,x,n,10);`):

elif nops([args])=1 and op(1,[args])=GuessFactor then
print(`GuessFactor(P,r): Given a polynomial P in r`):
print(`finds the polynomial f(x) such that the 2r-th moment of`):
print(`exp(-x^2/2)*f(x) is P. For example, try:`):
print(`GuessFactor(1,r,x);`):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(L,n): guesses a polynomial of degree d in n for`):
print(` the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: `):
print(`GuessPol([seq(i,i=1..10)],n);`):



elif nops([args])=1 and op(1,[args])=Moms then
print(`Moms(f,x,n,L): the first L moments (about the mean)`):
print(`of the distribution obtained by n throws of a die`):
print(`whose prob. generating function is f(x)`):
print(`For example, for the first 10 moments `):
print(`for tossing a loaded coin whose`):
print(` probability of Head is p, and the r.v. is [number of Heads], type: ` ):
print(`Moms(p*x+1-p,x,n,10);`):


elif nops([args])=1 and op(1,[args])=NorMomsE then
print(`NorMomsE(f,x,n,L): the first L normalized `):
print(`even moments (about the mean)`):
print(` i.e. the list whose r-th entry is the (2r)-th moment`):
print(`divided by the r-th power of the second moment`):
print(`divided by (2r)!/2^r/r!`):
print(`of the distribution obtained by n throws of a die`):
print(`whose prob. generating function is`):
print(`f(x). For example, try (for tossing a loaded coin whose`):
print(`probability of Head is p, and the r.v. is "number of Heads") type: `):
print(`NorMomsE(p*x+1-p,x,n,10);`):


elif nops([args])=1 and op(1,[args])=NorFaMomsE then
print(`NorFaMomsE(f,x,n,L): the first L normalized `):
print(`even Factorial moments (about the mean)`):
print(` i.e. the list whose r-th entry is the (2r)-th moment`):
print(`divided by the r-th power of the second moment`):
print(`divided by (2r)!/2^r/r!`):
print(`of the distribution obtained by n throws of a die`):
print(`whose prob. generating function is`):
print(`f(x). For example, try (for tossing a loaded coin whose`):
print(`probability of Head is p, and the r.v. is "number of Heads") type: `):
print(`NorFaMomsE(p*x+1-p,x,n,10);`):

elif nops([args])=1 and op(1,[args])=NorFaMomsO then
print(`NorFaMomsO(f,x,n,L): the first L normalized `):
print(`odd moments (about the mean)`):
print(` i.e. the list whose r-th entry is the (2r+1)-th Factorial moment`):
print(`divided by the r-th power of the second moment`):
print(`divided by (2r)!/2^r/r!`):
print(`of the distribution obtained by n throws of a die`):
print(`whose prob. generating function is`):
print(`f(x). For example, try (for tossing a loaded coin whose`):
print(`probability of Head is p, and the r.v. is "number of Heads") type: `):
print(`NorFaMomsO(p*x+1-p,x,n,10);`):

elif nops([args])=1 and op(1,[args])=NorMomsO then
print(`NorMomsO(f,x,n,L): the first L normalized `):
print(`odd moments (about the mean)`):
print(` i.e. the list whose r-th entry is the (2r+1)-th moment`):
print(`divided by the r-th power of the second moment`):
print(`divided by (2r)!/2^r/r!`):
print(`of the distribution obtained by n throws of a die`):
print(`whose prob. generating function is`):
print(`f(x). For example, try (for tossing a loaded coin whose`):
print(`probability of Head is p, and the r.v. is "number of Heads") type: `):
print(`NorMomsO(p*x+1-p,x,n,10);`):

elif nops([args])=1 and op(1,[args])=Prob then
print(`Prob(f,x,n,a,b): Given a die whose p.g.f. is f(x)`):
print(`and an integer n, the exact probability that the sum`):
print(`of the throws is between a and b inclusive.`):
print(`For example, for the prob. of throwing a `):
print(`fair coin 200 times and getting between 90 and 110`):
print(`Heads, type:`):
print(` Prob((1+x)/2,x,200,90,110);`):


elif nops([args])=1 and op(1,[args])=ProveAsyFaMom then
print(`ProveAsyFaMom(f,x,s): Rigorously proves the`):
print(`conjectured asymptotic expression given by both`):
print(`AsyFaMomE(f,x,n,r,s) and AsyFaMomO(f,x,n,r,s),`):
print(`i.e., the asymptotic expression of order s for all the`):
print(`factorial moments of the random variable "total number of points"`):
print(`upon throwing a die whose prob. generating function is`):
print(`f(x). For example, try:`):
print(`ProveAsyFaMom((1+x)/2,x,2);`):

elif nops([args])=1 and op(1,[args])=ProveAsyFaMomE then
print(`ProveAsyFaMomE(f,x,s): Rigorously proves the`):
print(`conjectured asymptotic expression given by AsyFaMomE(f,x,n,r,s),`):
print(`i.e., the asymptotic expression of order s for the even`):
print(`factorial moments of the random variable "total number of points"`):
print(`upon throwing a die whose prob. generating function is`):
print(`f(x). For example, try:`):
print(`ProveAsyFaMomE((1+x)/2,x,2);`):

elif nops([args])=1 and op(1,[args])=ProveAsyFaMomEverbose then
print(`ProveAsyFaMomEverbose(f,x,s): verbose form of`):
print(` ProveAsyFaMomE(f,x,s) `):
print(`For example, try:`):
print(`ProveAsyFaMomEverbose((1+x)/2,x,2);`):


elif nops([args])=1 and op(1,[args])=ProveAsyFaMomO then
print(`ProveAsyFaMomO(f,x,s): Rigorously proves the`):
print(`conjectured asymptotic expression given by AsyFaMomO(f,x,n,r,s),`):
print(`i.e., the asymptotic expression of order s for the odd`):
print(`factorial moments of the random variable "total number of points"`):
print(`upon throwing a die whose prob. generating function is`):
print(`f(x). For example, try:`):
print(`ProveAsyFaMomO((1+x)/2,x,2);`):

elif nops([args])=1 and op(1,[args])=ProveAsyFaMomOverbose then
print(`ProveAsyFaMomOverbose(f,x,s): verbose form of`):
print(` ProveAsyFaMomO(f,x,s)  `):
print(`For example, try:`):
print(`ProveAsyFaMomOverbose((1+x)/2,x,2);`):

elif nops([args])=1 and op(1,[args])=Stir2 then
print(`Stir2(r,i): The polynomial in r that is Stirling2(r,r-i).`):
print(`For example try: Stir2(r,2):`):


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

with(combinat):




##begin guessing polynomials
#GuessPol1(L,d,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i+1] for i=1..nops(L)
#For example, try: 
#GuessPol1([seq(i,i=1..10)],1,n);
GuessPol1:=proc(L,d,n) local P,i,a,eq,var:
if d>nops(L)-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=1..nops(L))}:

var:=solve(eq,var):

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

subs(var,P):

end:

#GuessPol(L,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i] for i=1..nops(L) for example, try: 
#GuessPol([seq(i,i=1..10)],n);
GuessPol:=proc(L,n) local d,gu:

for d from 0 to nops(L)-2 do
 gu:=GuessPol1(L,d,n):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:
##end guessing polynomials


#IsSymm(f,x): is the polynomial f symmetric?
IsSymm:=proc(f,x) local d1,d2,i:
d1:=ldegree(f,x):
d2:=degree(f,x):
if {seq(normal(coeff(f,x,d1+i)-coeff(f,x,d2-i)),i=0..d2-d1)}={0} then
 RETURN(true):
else
   RETURN(false):
fi:
end:





#Moms(f,x,n,L): the first L moments (about the mean)
#of the distribution obtained by n throws of a die
#whose prob. generating function is
#f(x). For example, try:
#Moms(p*x+1-p,x,n,10);
Moms:=proc(f,x,n,L) local gu,i,av,F:
option remember:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

gu:=[]:
av:=subs(x=1,diff(f,x)):
F:=(f/x^av)^n:

for i from 1 to L do
 F:=expand(x*diff(F,x)):
 gu:=[op(gu),normal(subs(x=1,F))]:
od:

gu:
end:

#NorMomsE(f,x,n,L): the first L normalized 
#even moments (about the mean)
#i.e. the list whose r-th entry is the (2r)-th moment
#divided by the r-th power of the second moment
#divided by (2r)!/2^r/r!
#of the distribution obtained by n throws of a die
#whose prob. generating function is
#f(x). For example, try (for tossing a loaded coin whose
#probability of Head is p, and the r.v. is "number of Heads") type: 
#NorMomsE(p*x+1-p,x,n,10)
NorMomsE:=proc(f,x,n,L) local gu,lu,i,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

gu:=Moms(f,x,n,2*L):
lu:=gu[2]:
gu:=[seq(gu[2*i],i=1..L)]:
gu:=expand([seq(normal(gu[i]/lu^i/(2*i)!*2^i*i!),i=1..L)]):

[seq(add(factor(coeff(gu[i],n,-j))/n^j,j=0..-ldegree(gu[i],n)),i=1..L)]:
end:


#NorMomsO(f,x,n,L): the first L normalized 
#odd moments (about the mean)
#i.e. the list whose r-th entry is the (2r+1)-th moment
#divided by the r-th power of the second moment
#divided by (2r)!/2^r/r!
#of the distribution obtained by n throws of a die
#whose prob. generating function is
#f(x). For example, try (for tossing a loaded coin whose
#probability of Head is p, and the r.v. is "number of Heads") type: 
#NorMomsO(p*x+1-p,x,n,10)
NorMomsO:=proc(f,x,n,L) local gu,lu,i,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

gu:=Moms(f,x,n,2*L+1):
lu:=gu[2]:
gu:=[seq(gu[2*i+1],i=1..L)]:
gu:=expand([seq(normal(gu[i]/lu^i/(2*i)!*2^i*i!),i=1..L)]):

[seq(add(factor(coeff(gu[i],n,-j))/n^j,j=0..-ldegree(gu[i],n)),i=1..L)]:
end:



#AsyNorMomE(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomE(1-p+p*x,x,n,r,3);
AsyNorMomE:=proc(f,x,n,r,s) local gu,L,lu,i,ku,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

L:=3*s+6:
lu:=NorMomsE(f,x,n,L):
gu:=0:

for i from 0 to s do
 ku:=[seq(coeff(lu[j],n,-i),j=1..L)]:
 ku:=GuessPol(ku,r):
  if ku=FAIL then
     RETURN(FAIL):
  fi:
  gu:=gu+factor(ku)/n^i:
od:

gu:

end:


#AsyNorMomO(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomO(1-p+p*x,x,n,r,3);
AsyNorMomO:=proc(f,x,n,r,s) local gu,L,lu,i,ku,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

L:=3*s+6:
lu:=NorMomsO(f,x,n,L):
gu:=0:

for i from 0 to s do
 ku:=[seq(coeff(lu[j],n,-i),j=1..L)]:
 ku:=GuessPol(ku,r):
  if ku=FAIL then
     RETURN(FAIL):
  fi:
  gu:=gu+factor(ku)/n^i:
od:

gu:

end:



#AsyMomE(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomE(1-p+p*x,x,n,r,3);
AsyMomE:=proc(f,x,n,r,s) local gu,lu:
lu:=factor(Moms(f,x,n,2)[2]/n):
gu:=AsyNorMomE(f,x,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:

#AsyMomO(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomO(1-p+p*x,x,n,r,3);
AsyMomO:=proc(f,x,n,r,s) local gu,lu:
lu:=factor(Moms(f,x,n,2)[2]/n):
gu:=AsyNorMomO(f,x,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:



#GuessFactor(P,r): Given a polynomial P in r
#finds the polynomial f(x) such that the 2r-th moment of
#exp(-x^2/2)*f(x) is P. For example, try:
#GuessFactor(1,r,x);
GuessFactor:=proc(P,r,x) local a,d,i,j,eq,var,gu:
d:=degree(P,r):
gu:=expand(add(a[i]*mul(2*r+2*j+1,j=0..i-1),i=0..d)-P):
var:={seq(a[i],i=0..d)}:
eq:={seq(coeff(gu,r,i),i=0..d)}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:
subs(var,add(a[i]*x^(2*i),i=0..d)):
end:



#CorrectedGaussian(f,x,n,s): The 
#Corrected Gaussian to the limiting distribution
#of a symmetric p.d. using the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x and n are symbolic, but s is a numeric non-neg. integer 
#For example, try: CorrectedGaussian((1+x)/2,x,n,3);
CorrectedGaussian:=proc(f,x,n,s) local r,gu,lu,i:
option remember:

if not IsSymm(f,x) then
 print(`Not symmetric, case not yet implemented`):
 RETURN(FAIL):
fi:

gu:=AsyNorMomE(f,x,n,r,s):

lu:=add(GuessFactor(coeff(gu,n,-i),r,x)/n^i,i=0..s):

exp(-x^2/2)*lu:

end:


#Prob(f,x,n,a,b): Given a die whose p.g.f. is f(x)
#and an integer n, the exact probability that the sum
#of the throws is between a and b inclusive.
#For example, for the prob. of throwing a 
#fair coin 200 times and getting between 90 and 110
#Heads, type:
# Prob((1+x)/2,x,200,90,110);
Prob:=proc(f,x,n,a,b) local lu,i:
lu:=expand(f^n):
add(coeff(lu,x,i),i=a..b):
end:


#AppxProb(f,x,n,a,b,s): Given a die whose p.g.f. is f(x)
#and an integer n, the approximate probability that the sum
#of the throws is between a and b inclusive.
#using a corrected Gaussian of order s.
#(s=0 corresponds to the usual  Normal Appx.)
#For example, for the appx. prob. of throwing a 
#fair coin 200 times and getting between 90 and 110
#Heads, using a corrected Gaussian of order 1, type:
# AppxProb((1+x)/2,x,200,90,110,1);
AppxProb:=proc(f,x,n1,a,b,s) local av, lu,a1,b1,PHI,nor1,n:

if not IsSymm(f,x) then
 print(`Not symmetric, case not yet implemented`):
 RETURN(FAIL):
fi:

av:=n1*evalf(subs(x=1,diff(f,x))):

lu:=evalf(subs(n=n1,sqrt(Moms(f,x,n,2)[2]))):

a1:=((a-1/2)-av)/lu:

b1:=((b+1/2)-av)/lu:

PHI:=subs(n=n1,CorrectedGaussian(f,x,n,s)):

nor1:=evalf(int(PHI,x=-infinity..infinity)):


evalf(int(PHI,x=a1..b1))/nor1:

end:




#Compare(f,x,n1,a,b,s): Compares the errors in the approximations
#by the corrected Gaussians to the uncorrected ones
#for exaple try:
#Compare((1+x)/2,x,100,40,60,3);
Compare:=proc(f,x,n1,a,b,s) local lu,s1,meduyak,taut:
meduyak:=evalf(Prob(f,x,n1,a,b)):

if not IsSymm(f,x) then
 print(`Not symmetric, case not yet implemented`):
 RETURN(FAIL):
fi:

for s1 from 0 to s do
lu[s1]:=AppxProb(f,x,n1,a,b,s1):
od:

taut:=abs(meduyak-lu[0]):

[seq(abs(lu[s1]-meduyak)/taut,s1=1..s)]:

end:



###Start Factorial Moments

#FactorialMoms(f,x,n,L): the first L Factorial  moments (about the mean)
#of the distribution obtained by n throws of a die
#whose prob. generating function is
#f(x). For example, try:
#FactorialMoms(p*x+1-p,x,n,10);
FactorialMoms:=proc(f,x,n,L) local gu,i,av,F:
option remember:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

gu:=[]:
av:=subs(x=1,diff(f,x)):
F:=(f/x^av)^n:

for i from 1 to L do
 F:=expand(diff(F,x)):
 gu:=[op(gu),normal(subs(x=1,F))]:
od:

gu:
end:



#NorFaMomsE(f,x,n,L): the first L normalized 
#even Factorial moments (about the mean)
#i.e. the list whose r-th entry is the (2r)-th Factorial moment
#divided by the r-th power of the second moment
#divided by (2r)!/2^r/r!
#of the distribution obtained by n throws of a die
#whose prob. generating function is
#f(x). For example, try (for tossing a loaded coin whose
#probability of Head is p, and the r.v. is "number of Heads") type: 
#NorFaMomsE(p*x+1-p,x,n,10)
NorFaMomsE:=proc(f,x,n,L) local gu,lu,i,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

gu:=FactorialMoms(f,x,n,2*L):
lu:=gu[2]:
gu:=[seq(gu[2*i],i=1..L)]:
gu:=expand([seq(normal(gu[i]/lu^i/(2*i)!*2^i*i!),i=1..L)]):

[seq(add(factor(coeff(gu[i],n,-j))/n^j,j=0..-ldegree(gu[i],n)),i=1..L)]:
end:


#NorFaMomsO(f,x,n,L): the first L normalized 
#odd Factorial moments (about the mean)
#i.e. the list whose r-th entry is the (2r+1)-th moment
#divided by the r-th power of the second moment
#divided by (2r)!/2^r/r!
#of the distribution obtained by n throws of a die
#whose prob. generating function is
#f(x). For example, try (for tossing a loaded coin whose
#probability of Head is p, and the r.v. is "number of Heads") type: 
#NorFaMomsO(p*x+1-p,x,n,10)
NorFaMomsO:=proc(f,x,n,L) local gu,lu,i,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

gu:=FactorialMoms(f,x,n,2*L+1):
lu:=gu[2]:
gu:=[seq(gu[2*i+1],i=1..L)]:
gu:=expand([seq(normal(gu[i]/lu^i/(2*i)!*2^i*i!),i=1..L)]):

[seq(add(factor(coeff(gu[i],n,-j))/n^j,j=0..-ldegree(gu[i],n)),i=1..L)]:
end:






#AsyNorFaMomO(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th Factorial moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorFaMomO(1-p+p*x,x,n,r,3);
AsyNorFaMomO:=proc(f,x,n,r,s) local gu,L,lu,i,ku,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

L:=3*s+6:
lu:=NorFaMomsO(f,x,n,L):
gu:=0:

for i from 0 to s do
 ku:=[seq(coeff(lu[j],n,-i),j=1..L)]:
 ku:=GuessPol(ku,r):
  if ku=FAIL then
     RETURN(FAIL):
  fi:
  gu:=gu+factor(ku)/n^i:
od:

gu:

end:



#AsyFaMomE(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized Factorial (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyFaMomE(1-p+p*x,x,n,r,3);
AsyFaMomE:=proc(f,x,n,r,s) local gu,lu:
lu:=factor(Moms(f,x,n,2)[2]/n):
gu:=AsyNorFaMomE(f,x,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:



#AsyNorFaMomE(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th Factorial moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorFaMomE(1-p+p*x,x,n,r,3);
AsyNorFaMomE:=proc(f,x,n,r,s) local gu,L,lu,i,ku,j:
if normal(subs(x=1,f))<>1 then
 print(f, `is not a prob. generating function`):
 RETURN(FAIL):
fi:

L:=3*s+6:
lu:=NorFaMomsE(f,x,n,L):
gu:=0:

for i from 0 to s do
 ku:=[seq(coeff(lu[j],n,-i),j=1..L)]:
 ku:=GuessPol(ku,r):
  if ku=FAIL then
     RETURN(FAIL):
  fi:
  gu:=gu+factor(ku)/n^i:
od:

gu:

end:

#AsyMomFaO(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized Factorial (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyFaMomO(1-p+p*x,x,n,r,3);
AsyFaMomO:=proc(f,x,n,r,s) local gu,lu:
lu:=factor(Moms(f,x,n,2)[2]/n):
gu:=AsyNorFaMomO(f,x,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:



#ProveAsyFaMomE(f,x,s): Rigorously proves the
#conjectured asymptotic expression given by AsyFaMomE(f,x,n,r,s),
#i.e., the asymptotic expression of order s for the even
#factorial moments of the random variable "total number of points"
#upon throwing a die whose prob. generating function is
#f(x). For example, try:
#ProveAsyFaMomE((1+x)/2,x,2);
ProveAsyFaMomE:=proc(f,x,s) local n,r,z,guE,guO,s1,lu,luE,
P,av,i,F,y,ka:

if normal(subs(x=1,f))<>1 then
  print(`Not a prob. gen. function`):
 RETURN(FAIL):
fi:

av:=normal(subs(x=1,diff(f,x))):
F:=f/x^av:

F:=subs(x=1+z,F):
F:=taylor(F,z=0,2*s+3):
P:=[seq(normal(coeff(F,z,i)),i=1..2*s+2)]:

lu:=factor(Moms(f,x,n,2)[2]/n):
guE:=AsyNorFaMomE(f,x,n,r,s):
guO:=AsyNorFaMomO(f,x,n,r,s):

guE:=subs({seq(1/n^s1=y^s1, s1=1..s+1)},guE):
guO:=subs({seq(1/n^s1=y^s1, s1=1..s+1)} ,guO):

luE:=(1+y)^r*subs(y=y/(1+y),guE)-guE:

for s1 from 1 to s do
luE:=luE-P[2*s1]*subs(r=r-s1,guE)/lu^s1*2^s1*simplify(r!/(r-s1)!)*y^s1:
od:


for s1 from 1 to s-1 do
luE:=luE-P[2*s1+1]*subs(r=r-s1-1,guO)/lu^(s1+1)*
2^(s1+1)*simplify(r!/(r-s1-1)!)*y^(s1+1)/(2*(r-s1-1)+1):
od:


luE:=taylor(luE,y=0,s+2):

ka:=[seq(normal(coeff(luE,y,i)),i=0..s)]:
if ka=[0$(s+1)] then
 RETURN(true):
else
  RETURN(ka):
fi:



end:



#ProveAsyFaMomO(f,x,s): Rigorously proves the
#conjectured asymptotic expression given by AsyFaMomO(f,x,n,r,s),
#i.e., the asymptotic expression of order s for the odd
#moments of the random variable "total number of points"
#upon throwing a die whose prob. generating function is
#f(x). For example, try:
#ProveAsyFaMomO((1+x)/2,x,2);
ProveAsyFaMomO:=proc(f,x,s) local n,r,z,guE,guO,s1,lu,luO,
P,av,i,F,y,ka:

if normal(subs(x=1,f))<>1 then
  print(`Not a prob. gen. function`):
 RETURN(FAIL):
fi:

av:=normal(subs(x=1,diff(f,x))):
F:=f/x^av:

F:=subs(x=1+z,F):
F:=taylor(F,z=0,2*s+3):
P:=[seq(normal(coeff(F,z,i)),i=1..2*s+2)]:

lu:=factor(Moms(f,x,n,2)[2]/n):
guE:=AsyNorFaMomE(f,x,n,r,s):
guO:=AsyNorFaMomO(f,x,n,r,s):

guE:=subs({seq(1/n^s1=y^s1, s1=1..s+1)},guE):
guO:=subs({seq(1/n^s1=y^s1, s1=1..s+1)} ,guO):

luO:=(1+y)^r*subs(y=y/(1+y),guO)-guO:

for s1 from 1 to s do
luO:=luO-P[2*s1]*subs(r=r-s1,guO)/lu^s1*2^s1*simplify(r!/(r-s1)!)*y^s1
*(2*r+1)/(2*(r-s1)+1):
od:


for s1 from 1 to s do
luO:=luO-P[2*s1+1]*subs(r=r-s1,guE)/lu^s1*
2^s1*simplify(r!/(r-s1)!)*y^s1*(2*r+1):
od:


luO:=taylor(luO,y=0,s+2):

ka:=[seq(normal(coeff(luO,y,i)),i=0..s)]:
if ka=[0$(s+1)] then
 RETURN(true):
else
  RETURN(ka):
fi:



end:

ProveAsyFaMom:=proc(f,x,s):
ProveAsyFaMomE(f,x,s) and ProveAsyFaMomO(f,x,s):
end:

#Stir2(r,i): The polynomial in r that is Stirling2(r,r-i)
#For example try: Stir2(r,2):
Stir2:=proc(r,i) local gu,j:
gu:=[seq(stirling2(j,j-i),j=1..3*i+2)]:
factor(GuessPol(gu,r)):
end:





#AsyNorMomEviaFa(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#via the Factorial moments
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomEviaFa(1-p+p*x,x,n,r,3);
AsyNorMomEviaFa:=proc(f,x,n,r,s) local gu,i,lu:
lu:=factor(Moms(f,x,n,2)[2]/n):
gu:=0:

for i from 0 to s do
gu:=gu+subs(r=2*r,Stir2(r,2*i))*subs(r=r-i,AsyNorFaMomE(f,x,n,r,s))
*simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i/lu^i/n^i:
:
od:

for i from 0 to s-1 do
gu:=gu+subs(r=2*r,Stir2(r,2*i+1))*subs(r=r-i-1,AsyNorFaMomO(f,x,n,r,s))
*simplify((2*r-2*i-2)!/(2*r)!)*simplify(r!/(r-i-1)!)*2^(i+1)/lu^(i+1)/
n^(i+1):
:
od:

add(factor(coeff(gu,n,-i))/n^i,i=0..s):

end:



#AsyNorMomOviaFa(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#via the Factorial moments
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomOviaFa(1-p+p*x,x,n,r,3);
AsyNorMomOviaFa:=proc(f,x,n,r,s) local gu,i,lu,vu:
lu:=factor(Moms(f,x,n,2)[2]/n):
gu:=0:



vu:=AsyNorFaMomO(f,x,n,r,s):
for i from 0 to s do
gu:=gu+subs(r=2*r+1,Stir2(r,2*i))*
subs(r=r-i,vu)
*simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i/lu^i/n^i:
:
od:


for i from 0 to s do
gu:=gu+subs(r=2*r+1,Stir2(r,2*i+1))*subs(r=r-i,AsyNorFaMomE(f,x,n,r,s))
*simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i/lu^i/n^i:
:
od:


add(factor(coeff(gu,n,-i))/n^i,i=0..s):

end:







#ProveAsyFaMomEverbose(f,x,s): Human readable
#statment and proof 
#for the asymptotic expression given by AsyFaMomE(f,x,n,r,s),
#i.e., the asymptotic expression of order s for the even
#factorial moments of the random variable "total number of points"
#upon throwing a die whose prob. generating function is
#f(x). For example, try:
#ProveAsyFaMomEverbose((1+x)/2,x,2);
ProveAsyFaMomEverbose:=proc(f,x,s) local n,r,z,guE,guO,s1,lu,luE,
P,av,i,F,y,ka,F1:

if normal(subs(x=1,f))<>1 then
  print(`Not a prob. gen. function`):
 RETURN(FAIL):
fi:

print(`Theorem:`):
print(`Consider a die whose probability generating function is`):
print(f):
print(`Consider the random variable: total score after n throws`):
print(`and Let A_r(n) be its r-th factorial moment,`):
print(` about the mean.`):
print(`The asymptotic expressions for the even, A_2r(n) and`):
print(`odd, A_(2r+1)(n) , are as follows.`):
print(A_(2*r)(n)=AsyFaMomE(f,x,n,r,s)):
print(`while the factorial odd, (2r+1)-th, moment is given by`):
print(A_(2*r+1)(n)=AsyFaMomO(f,x,n,r,s)):
print(`Proof`):
print(`The average outcome of each indiviual throw is the `):
print(`first derivative of `, f, `evaluated at x=1, which is`):
av:=normal(subs(x=1,diff(f,x))):
print(av):
F1:=f/x^av:
print(`Dividing by `, x^av, ` gives that that the centralized p.g.f is`):
print(`F:=`, F(x)=F1):
print(`and the p.g.f. for n throws is`):
print(F(x)^n=F1^n):
print(`Recall that the factorial moments satisfy`):
print(Sum(A_r(n)*z^r/r!,r=0..infinity)=F(1+z)^n):
print(`Hence we have`):
print(Sum(A_r(n+1)*z^r/r!,r=0..infinity)=Sum(A_r(n)*z^r/r!,r=0..infinity)*
F(1+z)):
F1:=subs(x=1+z,F1):
F1:=taylor(F1,z=0,2*s+3):
print(`Let's expand , F(1+z), as a power series around z=0, up`):
print(`to order`, 2*s+2, `getting:`):
print(F1):
P:=[seq(normal(coeff(F1,z,i)),i=1..2*s+2)]:
print(`The list of first`, 2*s+2, `coefficients, let's call them  C, is:`):
print(P):
print(`Comparing coefficients in the above recurrence`):
print(`we have, to order`, s, `the recurrence`):
print(A_r(n+1)-A_r(n)=add(C[i]*A_(r-i)(n),i=1..2*s)):
print(`which spells out to:`):
print(A_r(n+1)-A_r(n)=add(P[i]*A_(r-i)(n),i=1..2*s)):

lu:=factor(Moms(f,x,n,2)[2]/n):

print(`Spelling it out separately for the even and odd cases, we have`):
print(
A_(2*r)(n+1)-A_(2*r)(n)=add(P[2*i]*A_(2*r-2*i)(n),i=1..s)+add(P[2*i+1]*A_(2*r-2*i-1)(n),i=1..s-1)):
print():
print(
A_(2*r+1)(n+1)-A_(2*r+1)(n)=add(P[2*i]*A_(2*r-2*i+1)(n),i=1..s)+add(P[2*i+1]*A_(2*r-2*i)(n),i=1..s)):

print(`Let's write, for the even case `):
print(A_(2*r)(n)=(2*r)!/r!/2^r*lu^n*n^r*B_r(n)):
print(`and let's write`):
print(A_(2*r+1)(n)=(2*r)!/r!/2^r*lu^n*n^r*C_r(n)):
print(`for the odd case`):

print(`We have to prove that` ):
print(B_r(n)=AsyNorFaMomE(f,x,n,r,s)+O(1/n^(s+1))):
print(C_r(n)=AsyNorFaMomO(f,x,n,r,s)+O(1/n^(s+1))):


print(`We have to prove that B_r(n), C_r(n) satisfy`):
print(`up to order`, 1/n^s):
print(
((n+1)/n)^r*
B_r(n+1)-B_r(n)=add(P[2*i]*B_(r-i)(n)/lu^i*2^i*simplify(r!/(r-i)!),i=1..s)
+add(P[2*i+1]/n^i*C_(r-i-1)(n)
/lu^(i+1)*2^(i+1)*simplify(r!/(r-i-1)!)/n^(i+1)
,i=1..s-1)
):
print():
print(
((n+1)/n)^r*C_r(n+1)-C_(r)(n)=add(P[2*i]*C_(r-i)(n)*
lu^i*2^i*simplify(r!/(r-i)!)*(2*r+1)/(2*(r-i)+1)/n^i
,i=1..s)
+add(P[2*i+1]*B_(r-i)(n)/lu^i*2^i*simplify(r!/(r-i)!)
,i=1..s)*(2*r+1)):


guE:=AsyNorFaMomE(f,x,n,r,s):
guO:=AsyNorFaMomO(f,x,n,r,s):

guE:=subs({seq(1/n^s1=y^s1, s1=1..s+1)},guE):
guO:=subs({seq(1/n^s1=y^s1, s1=1..s+1)} ,guO):

luE:=(1+y)^r*subs(y=y/(1+y),guE)-guE:

for s1 from 1 to s do
luE:=luE-P[2*s1]*subs(r=r-s1,guE)/lu^s1*2^s1*simplify(r!/(r-s1)!)*y^s1:
od:


for s1 from 1 to s-1 do
luE:=luE-P[2*s1+1]*subs(r=r-s1-1,guO)/lu^(s1+1)*
2^(s1+1)*simplify(r!/(r-s1-1)!)*y^(s1+1)/(2*(r-s1-1)+1):
od:


luE:=taylor(luE,y=0,s+2):

print(`Plugging-everything in, and expanding up to`, 1/n^s):
print(`we get that everything is indeed`):

ka:=[seq(normal(coeff(luE,y,i)),i=0..s)]:
if ka=[0$(s+1)] then
print(`true, QED`):
else
 print(`false, too bad!`):
fi:



end:



#################New stuff




#FactorialMomsG(P,n,L): the first L Factorial  moments (about the mean)
#of the distribution obtained by n throws of a die
#whose p.d.f. has factorial moments P[1]=0, P[2], P[3], ...  .
#For example, try:
#FactorialMomsG(P,n,10);
FactorialMomsG:=proc(P,n,L) local i,av,F,z:
option remember:

F:=(1+add(P[i]/i!*z^i,i=2..L))^n:
F:=taylor(F,z=0,L+1):
[seq(expand(i!*coeff(F,z,i)),i=1..L)]:

end:


#NorFaMomsEg(P,x,n,L): the first L normalized 
#even Factorial moments (about the mean)
#i.e. the list whose r-th entry is the (2r)-th Factorial moment
#divided by the r-th power of the second moment
#divided by (2r)!/2^r/r!
#of the distribution obtained by n throws of a 
#proverbial die
#whose factorial moments are P[2],P[3] (given symbolically)
#For example, try 
#NorFaMomsEg(P,x,n,10)
NorFaMomsEg:=proc(P,n,L) local gu,lu,i,j:

gu:=FactorialMomsG(P,n,2*L):
lu:=gu[2]:
gu:=[seq(gu[2*i],i=1..L)]:
gu:=expand([seq(normal(gu[i]/lu^i/(2*i)!*2^i*i!),i=1..L)]):

[seq(add(factor(coeff(gu[i],n,-j))/n^j,j=0..-ldegree(gu[i],n)),i=1..L)]:
end:

#AsyNorFaMomEg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th Factorial moment
#of the r.v. obtained by repeating n times throwing
# a proverbial die (with mean 0), whose factorial moments are P[2],P[3]
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorFaMomEg(P,n,r,3);
AsyNorFaMomEg:=proc(P,n,r,s) local gu,L,lu,i,ku,j:
option remember:
L:=3*s+6:
lu:=NorFaMomsEg(P,n,L):
gu:=0:

for i from 0 to s do
 ku:=[seq(coeff(lu[j],n,-i),j=1..L)]:
 ku:=GuessPol(ku,r):
  if ku=FAIL then
     RETURN(FAIL):
  fi:
  gu:=gu+factor(ku)/n^i:
od:

gu:

end:


#NorFaMomsOg(P,x,n,L): the first L normalized 
#even Factorial moments (about the mean)
#i.e. the list whose r-th entry is the (2r+1)-th Factorial moment
#divided by the r-th power of the second moment
#divided by (2r)!/2^r/r!
#of the distribution obtained by n throws of a 
#proverbial die
#whose factorial moments are P[2],P[3] (given symbolically)
#For example, try 
#NorFaMomsOg(P,x,n,10)
NorFaMomsOg:=proc(P,n,L) local gu,lu,i,j:
option remember:
gu:=FactorialMomsG(P,n,2*L+1):
lu:=gu[2]:
gu:=[seq(gu[2*i+1],i=1..L)]:
gu:=expand([seq(normal(gu[i]/lu^i/(2*i)!*2^i*i!),i=1..L)]):

[seq(add(factor(coeff(gu[i],n,-j))/n^j,j=0..-ldegree(gu[i],n)),i=1..L)]:
end:

#AsyNorFaMomOg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th Factorial moment
#of the r.v. obtained by repeating n times throwing
# a proverbial die (with mean 0), whose factorial moments are P[2],P[3]
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorFaMomOg(P,n,r,3);
AsyNorFaMomOg:=proc(P,n,r,s) local gu,L,lu,i,ku,j:
option remember:
L:=3*s+6:
lu:=NorFaMomsOg(P,n,L):
gu:=0:

for i from 0 to s do
 ku:=[seq(coeff(lu[j],n,-i),j=1..L)]:
 ku:=GuessPol(ku,r):
  if ku=FAIL then
     RETURN(FAIL):
  fi:
  gu:=gu+factor(ku)/n^i:
od:

gu:

end:


#PgfToF(f,x,L): the first L factorial moments
#(about the mean)
#of a discrete p.d. whose pgf is f(x)
#For example, try:
#PgfToF(p*x+1-p,x,10);
PgfToF:=proc(f,x,L) local i,z,F:
if normal(subs(x=1,f))<>1 then
 print(`Not a pgf`):
 RETURN(FAIL):
fi:
F:=normal(f/x^subs(x=1,diff(f,x))):
F:=subs(x=1+z,F):
F:=taylor(F,z=0,L+1):
[seq(normal(i!*coeff(F,z,i)),i=1..L)]:
end:





#AsyNorFaMomO1(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th Factorial moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorFaMomO1(1-p+p*x,x,n,r,3);
AsyNorFaMomO1:=proc(f,x,n,r,s) local gu,mu,P,i1:

gu:=AsyNorFaMomOg(P,n,r,s):
mu:=PgfToF(f,x,3*s+6):
gu:=subs({seq(P[i1]=mu[i1],i1=1..3*s+6)},gu):
add(factor(coeff(gu,n,-i1))/n^i1,i1=0..s):
end:



#AsyNorFaMomE1(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th Factorial moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorFaMomE1(1-p+p*x,x,n,r,3);
AsyNorFaMomE1:=proc(f,x,n,r,s) local gu,mu,P,i1:

gu:=AsyNorFaMomEg(P,n,r,s):
mu:=PgfToF(f,x,3*s+6):
gu:=subs({seq(P[i1]=mu[i1],i1=1..3*s+6)},gu):

add(factor(coeff(gu,n,-i1))/n^i1,i1=0..s):
end:


#AsyFaMomEg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized Factorial (2r)-th moment
#of the r.v. obtained by repeating n times a p.d. 
#whose factorial moments are P[2],P[3],
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyFaMomEg(P,n,r,3);
AsyFaMomEg:=proc(P,n,r,s) local gu,lu:
lu:=P[2]:
gu:=AsyNorFaMomEg(P,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:

#AsyFaMomOg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized Factorial (2r+1)-th moment
#of the r.v. obtained by repeating n times a p.d. 
#whose factorial moments are P[2],P[3],
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyFaMomOg(P,n,r,3);
AsyFaMomOg:=proc(P,n,r,s) local gu,lu:
lu:=P[2]:
gu:=AsyNorFaMomOg(P,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:




#AsyNorMomEg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#via the Factorial moments
#For example, try: AsyNorMomEg(P,n,r,3);
AsyNorMomEg:=proc(P,n,r,s) local gu,i,lu:
lu:=P[2]:
gu:=0:

for i from 0 to s do
gu:=gu+subs(r=2*r,Stir2(r,2*i))*subs(r=r-i,AsyNorFaMomEg(P,n,r,s))
*simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i/lu^i/n^i:
:
od:

for i from 0 to s-1 do
gu:=gu+subs(r=2*r,Stir2(r,2*i+1))*subs(r=r-i-1,AsyNorFaMomOg(P,n,r,s))
*simplify((2*r-2*i-2)!/(2*r)!)*simplify(r!/(r-i-1)!)*2^(i+1)/lu^(i+1)/
n^(i+1):
:
od:

add(factor(coeff(gu,n,-i))/n^i,i=0..s):

end:



#AsyNorMomOg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#via the Factorial moments
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomOg(P,n,r,3);
AsyNorMomOg:=proc(P,n,r,s) local gu,i,lu,vu:
lu:=P[2]:
gu:=0:



vu:=AsyNorFaMomOg(P,n,r,s):
for i from 0 to s do
gu:=gu+subs(r=2*r+1,Stir2(r,2*i))*
subs(r=r-i,vu)
*simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i/lu^i/n^i:
:
od:


for i from 0 to s do
gu:=gu+subs(r=2*r+1,Stir2(r,2*i+1))*subs(r=r-i,AsyNorFaMomEg(P,n,r,s))
*simplify((2*r-2*i)!/(2*r)!)*simplify(r!/(r-i)!)*2^i/lu^i/n^i:
:
od:


add(factor(coeff(gu,n,-i))/n^i,i=0..s):

end:


#AsyMomEg(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a proverbial die whose factorial moments are given by P[2],P[3] etc.
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomEg(P,n,r,3);
AsyMomEg:=proc(P,n,r,s) local gu,lu:
lu:=P[2]:
gu:=AsyNorMomEg(P,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:

#AsyMomO(P,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a proverbial die whose factorial moments are P[2],P[3], ...
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomOg(P,n,r,3);
AsyMomOg:=proc(P,n,r,s) local gu,lu:
lu:=P[2]:
gu:=AsyNorMomOg(P,n,r,s):
n^r*lu^r*(2*r)!/r!/2^r*gu:
end:



#AsyMomE1(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the  (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomE1(1-p+p*x,x,n,r,3);
AsyMomE1:=proc(f,x,n,r,s) local gu,i1,P,mu:
gu:=AsyMomEg(P,n,r,s):
mu:=PgfToF(f,x,3*s+6):
gu:=subs({seq(P[i1]=mu[i1],i1=1..3*s+6)},gu):
end:

#AsyMomO1(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the  (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomO1(1-p+p*x,x,n,r,3);
AsyMomO1:=proc(f,x,n,r,s) local gu,i1,P,mu:
gu:=AsyMomOg(P,n,r,s):
mu:=PgfToF(f,x,3*s+6):
gu:=subs({seq(P[i1]=mu[i1],i1=1..3*s+6)},gu):
end:



#AsyNorMomE1(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomE1(1-p+p*x,x,n,r,3);
AsyNorMomE1:=proc(f,x,n,r,s) local gu,P,mu,i,i1:
gu:=AsyNorMomEg(P,n,r,s):
mu:=PgfToF(f,x,3*s+6):
gu:=subs({seq(P[i1]=mu[i1],i1=1..3*s+6)},gu):
add(factor(coeff(gu,n,-i))/n^i,i=0..s):
end:

#AsyNorMomO1(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomO1(1-p+p*x,x,n,r,3);
AsyNorMomO1:=proc(f,x,n,r,s) local gu,P,mu,i,i1:
gu:=AsyNorMomOg(P,n,r,s):
mu:=PgfToF(f,x,3*s+6):
gu:=subs({seq(P[i1]=mu[i1],i1=1..3*s+6)},gu):
add(factor(coeff(gu,n,-i))/n^i,i=0..s):
end:






#FactorialInTermsOfM(M,L): Expressing the first L factorial
#moments in terms of the moments M[1],M[2], ...
#For example, try: 
#FactorialInTermsOfM(M,5);
FactorialInTermsOfM:=proc(M,L) local r,gu,i,lu,x:
lu:=x:
gu:=[0]:

for r from 2 to L do
 lu:=expand(lu*(x-r+1)):
  gu:=[op(gu),add(coeff(lu,x,i)*M[i],i=2..degree(lu,x))]:
od:
gu:
end:

 







#AsyNorMomEgM(M,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a 
#via the Factorial moments
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomEg(P,n,r,3);
AsyNorMomEgM:=proc(M,n,r,s) local gu,P,vu:
option remember:
gu:=AsyNorMomEg(P,n,r,s):
vu:=FactorialInTermsOfM(M,2*s+4):
gu:=subs({seq(P[i]=vu[i],i=1..2*s+4)},gu):
add(factor(simplify(coeff(gu,n,-i)))/n^i,i=0..s):
end:


#AsyNorMomOgM(M,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a 
#via the Factorial moments
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomOg(P,n,r,3);
AsyNorMomOgM:=proc(M,n,r,s) local gu,P,vu:
option remember:
gu:=AsyNorMomOg(P,n,r,s):
vu:=FactorialInTermsOfM(M,2*s+4):
gu:=subs({seq(P[i]=vu[i],i=1..2*s+4)},gu):
add(factor(coeff(gu,n,-i))/n^i,i=0..s):
end:


#AsyMomEgM(M,n,r,s): The first s+1 terms in the
#asymptotic expansion of the  (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a 
#via the Factorial moments
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomEg(P,n,r,3);
AsyMomEgM:=proc(M,n,r,s) local gu,P,vu:

AsyNorMomEgM(M,n,r,s)*M[2]^r*n^r*(2*r)!/2^r/r!:
end:


#AsyMomOgM(M,n,r,s): The first s+1 terms in the
#asymptotic expansion of the  (2r+1)-th moment
#of the r.v. obtained by repeating n times throwing
# a 
#via the Factorial moments
#Here n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomOg(P,n,r,3);
AsyMomOgM:=proc(M,n,r,s) local gu,P,vu:
AsyNorMomOgM(M,n,r,s)*M[2]^r*n^r*(2*r)!/2^r/r!:
end:











#AsyNorMomEd(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomE(1-p+p*x,x,n,r,3);
AsyNorMomEd:=proc(f,x,n,r,s) local gu,M,lu,i:
gu:=AsyNorMomEgM(M,n,r,s):
lu:=subs(n=1,Moms(f,x,n,4*s)):
gu:=subs({seq(M[i]=lu[i],i=1..4*s)},gu):
add(factor(coeff(gu,n,-i))/n^i,i=0..s):
end:

#AsyNorMomOd(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the normalized (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomE(1-p+p*x,x,n,r,3);
AsyNorMomOd:=proc(f,x,n,r,s) local gu,M,lu,i:
gu:=AsyNorMomOgM(M,n,r,s):

lu:=subs(n=1,Moms(f,x,n,4*s+4)):
gu:=subs({seq(M[i]=lu[i],i=1..4*s+4)},gu):
add(factor(coeff(gu,n,-i))/n^i,i=0..s):
end:


#AsyMomEd(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyMomEd(1-p+p*x,x,n,r,3);
AsyMomEd:=proc(f,x,n,r,s) local gu,M,lu,i:
gu:=AsyMomEgM(M,n,r,s):
lu:=subs(n=1,Moms(f,x,n,4*s)):
gu:=subs({seq(M[i]=lu[i],i=1..4*s)},gu):

end:

#AsyMomOd(f,x,n,r,s): The first s+1 terms in the
#asymptotic expansion of the (2r)-th moment
#of the r.v. obtained by repeating n times throwing
# a die whose probability generating function is f(x)
#Here x,n,r are symbolic, but s is a numeric non-neg. integer 
#For example, try: AsyNorMomE(1-p+p*x,x,n,r,3);
AsyMomOd:=proc(f,x,n,r,s) local gu,M,lu,i:
gu:=AsyMomOgM(M,n,r,s):

lu:=subs(n=1,Moms(f,x,n,4*s+4)):
gu:=subs({seq(M[i]=lu[i],i=1..4*s+4)},gu):

end:







####New stuff



#ProveAsyFaMomEg(s): Rigorously proves the
#conjectured asymptotic expression given by AsyFaMomEg(P,n,r,s),
#i.e., the asymptotic expression of order s for the even
#factorial moments of the random variable obtained
#by repeating n times the elementary r.v. 
#For example, try:
#ProveAsyFaMomEg(2);
ProveAsyFaMomEg:=proc(s) local n,r,z,guE,guO,s1,lu,luE,
P,av,i,F,y,ka:

lu:=P[2]:
guE:=AsyNorFaMomEg(P,n,r,s):
guO:=AsyNorFaMomOg(P,n,r,s):

guE:=subs({seq(1/n^s1=y^s1, s1=1..s+1)},guE):
guO:=subs({seq(1/n^s1=y^s1, s1=1..s+1)} ,guO):

luE:=(1+y)^r*subs(y=y/(1+y),guE)-guE:

for s1 from 1 to s do
luE:=luE-(P[2*s1]/(2*s1)!)*subs(r=r-s1,guE)/lu^s1*2^s1*
simplify(r!/(r-s1)!)*y^s1:
od:


for s1 from 1 to s-1 do
luE:=luE-(P[2*s1+1]/(2*s1+1)!)*subs(r=r-s1-1,guO)/lu^(s1+1)*
2^(s1+1)*simplify(r!/(r-s1-1)!)*y^(s1+1)/(2*(r-s1-1)+1):
od:


luE:=taylor(luE,y=0,s+2):

ka:=[seq(normal(coeff(luE,y,i)),i=0..s)]:
if ka=[0$(s+1)] then
 RETURN(true):
else
  RETURN(ka):
fi:



end:


#ProveAsyFaMomOg(s): Rigorously proves the
#conjectured asymptotic expression given by AsyFaMomOg(P,n,r,s),
#For example, try:
#ProveAsyFaMomOg(2);
ProveAsyFaMomOg:=proc(s) local n,r,z,guE,guO,s1,lu,luO,
P,av,i,F,y,ka:

lu:=P[2]:
guE:=AsyNorFaMomEg(P,n,r,s):
guO:=AsyNorFaMomOg(P,n,r,s):

guE:=subs({seq(1/n^s1=y^s1, s1=1..s+1)},guE):
guO:=subs({seq(1/n^s1=y^s1, s1=1..s+1)} ,guO):

luO:=(1+y)^r*subs(y=y/(1+y),guO)-guO:

for s1 from 1 to s do
luO:=luO-(P[2*s1]/(2*s1)!)*
subs(r=r-s1,guO)/lu^s1*2^s1*simplify(r!/(r-s1)!)*y^s1
*(2*r+1)/(2*(r-s1)+1):
od:


for s1 from 1 to s do
luO:=luO-(P[2*s1+1]/(2*s1+1)!)*subs(r=r-s1,guE)/lu^s1*
2^s1*simplify(r!/(r-s1)!)*y^s1*(2*r+1):
od:


luO:=taylor(luO,y=0,s+2):

ka:=[seq(normal(coeff(luO,y,i)),i=0..s)]:
if ka=[0$(s+1)] then
 RETURN(true):
else
  RETURN(ka):
fi:



end:

ProveAsyFaMomG:=proc(s):
ProveAsyFaMomEg(s) and ProveAsyFaMomOg(s):
end:





#SipurG(M,n,r,s): Proves a refined version of the
#Central Limit Theorem, complete with asymptotics
#and the fact that if the first few moments
#coincide with the momenets of the Normal Distribution
#then the convergence to the Normal Distribution is
#faster; Do:
#SipurG(M,n,r,2);
SipurG:=proc(M,n,r,s) local gu,i,lu,mu,gu1:
gu:=AsyNorMomEgM(M,n,r,s):
gu1:=gu:
print(`A generalization of the Central Limit Theorem`):
print(``):
print(`by Shalosh B. Ekhad `):
print(``):
print(`Consider the centralized  distribution`):
print(`obtained by summing n identically-distributed `):
print(`independent (abbreviated :i.i.d.) `):
print(`r.v. with finite moments, M[1]=0, M[2],M[3], ...`):
print(`The 2r-th moment, divided by the (2r)-th moment of`):
print(`of the Standard Normal Distribution, (2r)!/(r!2^r)`):
print(`divided by M[2]^r has`):
print(`has asymptotics up to  order`, 1/n^s):
print(`as follows. `):

print(gu):
print(`If that distribution is symmetric, and the second moment`):
print(`is normalized to 1, then it equals`):
gu:=subs({seq(M[2*i+1]=0,i=0..s),M[2]=1},gu):
print(gu):

for i from 1 to 1 do
print(`if we want the`, 1/n^i , ` term to vanish `):
mu:=coeff(gu,n,-i):
lu:=solve(mu,M[2*i+2]):
print(`we need the`, 2*i+2, ` -th moment to be`, lu):
gu:=subs(M[2*i+2]=lu,gu):
od:

for i from 2 to s do
print(`if we also want the`, 1/n^i, ` term to vanish `):
mu:=coeff(gu,n,-i):
lu:=solve(mu,M[2*i+2]):
lu:=simplify(lu):
print(`we need the`, 2*i+2, ` -th moment to be`, lu):
gu:=subs(M[2*i+2]=lu,gu):
od:

gu1:

end: