######################################################################
##BiVariateMoms.txt: Save this file as  BiVariateMoms.txt            #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read BiVariateMoms.txt<Enter>                      #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec.  24, 2016
 
print(`Created: Dec. 15,  2016`):
print(` This is BiVariateMoms.txt `):
print(`A  Maple package to generate explicit formulas for moments of combinatorial random variables whose bi-variate generating`):
print(`functions are rational `):
print(``):
print(`It one of the packages that accompany 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 checking procedures type ezraC();, 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):




ezraC:=proc()

if args=NULL then
 print(` The CHECKING procedures are`):
 print(`  CheckAvMomsAM `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AsySize, AsySizeT, AveAndMoms, AvMoms, AvMomsAM,  AvMomsAstraight,AvMomsAstraightE, FindOpe, FindOpe1,`):
 print(`  MDktw,  Pashet, PW, PWs, PWsLL, RatToAsy, SeqT, ShoreshG, TestLLL  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Alpha, AlphaAsy, AlphaFl, AppxAnk, AsyData, AsyDataFl, AvMomsA, AvMomsAM, AvMomsAfl, AvMomsAE, ExactAnk,fmkNdata  `):
 print(`  LogCurveN, Story  `):
 print(` `):

elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(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,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,Variance,3rdStandarizedMoment, ..., JthStandarized]]] `):
print(`where L0 is the asymptotic formula for the actual size of the n-th set, in terms of the`):
print(`growth constant, denoted by  beta, Av, Variance, Skewness^2, Kurtosis..., JthAlphaCoeffi. are the asymptotic expressions`):
print(`in terms of beta  and n for the average, variance ..., all the way to the stadarized J-th moment, BUT `):
print(`Note that the odd-standarized moments are squared.`):
print(`(also about the mean), followed by the minimal polynomial in t, such that beta is its largest root.`):
print(`Try:`):
print(`Alpha(1/(1-t-t^2*w),t,w,4,n,beta);`):


elif nops([args])=1 and op(1,[args])=AlphaAsy then
print(`AlphaAsy(F,t,w,J,n,beta,r): Like AlphaAsy(F,t,w,J,n,beta) (q.v.) but with an asymptotic expression in n to order r`):
print(`of the standarized 3rd- to J-th moments (for the odd ones its square).`):
print(` Try: `):
print(` AlphaAsy(1/(1-t-t^2*w),t,w,6,n,beta,2); `):

elif nops([args])=1 and op(1,[args])=AlphaFl then
print(`AlphaFl(F,t,w,J,n): Like Alpha(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 STANDARIZED moment`):
print(`Try:`):
print(`AlphaFl(1/(1-t-t^2*w),t,w,4,n);`):


elif nops([args])=1 and op(1,[args])=AppxAnk then
print(`AppxAnk(F,t,w,n1,k1): the approximation of the coeff. of t^n1*w^k1 in the Taylor expansion of the rational function F in t,w`):
print(`using asymptotic normality (assuming local limit rule).`). 
print(`It returns a pair. The first is the appx. using the normal distibution`):
print(`and the second is how many standard-dviations it is from the mean. If this is very large or small the approximation is not good.`):
print(`Try:`):
print(`AppxAnk(1/(1-t-t^2*w),t,w,100,50);`):

elif nops([args])=1 and op(1,[args])=AsyData then
print(`AsyData(F,t,w,beta): inputs a rational function F in the variables t and w such that the coeff. of t^n is the generating function`):
print(`according to some statistic of the n-th member of an sequence of combinatorial objects, outputs the triple`):
print(`[A,c1,c2,beta,pol] where c1 and c2 are expressions in terms of the algebraric number beta, (the largest positive root`):
print(`of the polynomial in beta, and the meaning of c1 and c2 are that the asymptotics of the average of the statistic is c1*n and `):
print(` the asymptotics of the variance is c2*n. The asymptotics for the cardinality of the n-th family is A*beta^n.  Try:`):
print(`AsyData(1/(1-t-w*t^2),t,w,beta);`):


elif nops([args])=1 and op(1,[args])=AsyDataFl then
print(`AsyDataFl(F,t,w): Like AsyDataFl(F,t,w,beta) (q.v.) but in floating-point. `):
print(`It returns the four numbers [A,beta,c1,c2] such that the cardinality of the n-th member in the sequence is `):
print(`asympototically, A*beta^n, the expectation of the statistic is asymptotically c1*n, and the variance is asymptotically c2*n.`):
print(`Try:`):
print(`AsyDataFl(1/(1-t-w*t^2),t,w);`):

elif nops([args])=1 and op(1,[args])=AsySize then
print(`AsySize(R,t,n,beta): inputs a rational function R of the variable t, that is the  generating function  of some sequence of positive`):
print(`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`):
print(`that is also returned (the denominator of R with t->1/t). Try:`):
print(`AsySize(1/(1-t-t^2),t,n,beta); `):


elif nops([args])=1 and op(1,[args])=AsySizeT then
print(`AsySizeT(R,t,beta): inputs a rational function R of the variable t, that is the  generating function  of some sequence of positive`):
print(`integers a(n). Returns [A,beta,PolSatisfyingbeta] where the `):
print(` asympototic formula for a(n) in the form A*beta^n, where beta is the largest root of a polynomial in t`):
print(`that is also returned (the denominator of R with t->1/t). The output is the triple [A,beta,POLinBeta].  Try:`):
print(`AsySizeT(1/(1-t-t^2),t,beta);`):

elif nops([args])=1 and op(1,[args])=AveAndMoms then
print(`AveAndMoms(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list whose first entry is the average `):
print(`(alias expectation, alias mean)`):
print(`followed by the variance, the third moment (about the mean) and`):
print(`so on, until the N-th moment (about the mean).`):
print(`If f(1) is not 1, than it first normalizes it by dividing`):
print(`by f(1) (if it is not 0) .`):
print(`For example, try:`):
print(`AveAndMoms(((1+x)/2)^100,x,4);`):

elif nops([args])=1 and op(1,[args])=AvMoms then
print(`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,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs positive integers n0 and J. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]]  (straight moments)`):
print(`Try:`):
print(`AvMoms(1/(1-t*(1+w)),t,w,10,2);`):
print(`AvMoms(1/(1-t-w*t^2),t,w,30,2);`):
print(`AvMoms(MDktw(2,t,w),t,w,10,2);`):

elif nops([args])=1 and op(1,[args])=AvMomsA then
print(`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,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]] `):
print(`where L0 is the asymptotic formula for the actual size of the n-th set, in terms of the`):
print(`growth constant, denoted by  beta, Av, 2ndMom, ..., JthMom are the asymptotic expressions`):
print(`in terms of beta  and n for the average, 2ndMom (ABOUT THE MEAN), ..., all the way to the J-th moment`):
print(`(also ABOUT THE MEAN), followed by the minimal polynomial in t, such that beta is its largest root.`):
print(`Try: `):
print(`AvMomsA(1/(1-t-t^2*w),t,w,4,n,beta);`):


elif nops([args])=1 and op(1,[args])=AvMomsAE then
print(`AvMomsAE(F,t,w,J,n,a): inputs a rational function F of variables t, and w, such that the coefficient of t^n is the weight-enumerator, in w,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]] `):
print(`where  the actual size is denoted by a(n), and Av,2ndMom, ..., JthMom`):
print(`are the expressions, in terms of n and a(n), a(n+1), etc.`):
print(`for the NUMERATORS of the expressions for the average, and the moments ABOUT THE MEAN, ..., all the way to the J-th moment`):
print(` (also  ABOUT THE MEAN). To get the actual expression for the j-th moment, divide by a(n)^j. The second output is  the generating function of a(n) `):
print(`with respect to  the variable`,  t):
print(` Try: `):
print(`AvMomsAE(1/(1-t-t^2*w),t,w,4,n,a);`):

elif nops([args])=1 and op(1,[args])=AvMomsAM then
print(`AvMomsAM(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,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs positive integers n0 and J. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]]  with the sequence of numbers`):
print(`for the moments about the mean`):
print(`AvMomsAM(1/(1-t*(1+w)),t,w,10,2);`):
print(`AvMomsAM(MDktw(2,t,w),t,w,10,2);`):


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])=AvMomsAstraight then
print(`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,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]] `):
print(`where L0 is the asymptotic formula for the actual size of the n-th set, in terms of the`):
print(`growth constant, denoted by  beta, Av, 2ndMom, ..., JthMom are the asymptotic expressions`):
print(`in terms of beta  and n for the average, 2ndMom (NOT about the mean), ..., all the way to the J-th moment`):
print(`(also not about the mean), followed by the minimal polynomial in t, such that beta is its largest root.`):
print(`Try: `):
print(`AvMomsAstraight(1/(1-t-t^2*w),t,w,4,n,beta);`):



elif nops([args])=1 and op(1,[args])=AvMomsAstraightE then
print(`AvMomsAstraightE(F,t,w,J,n,a): inputs a rational function F of variables t, and w, such that the coefficient of t^n is the weight-enumerator, in w,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs a positive integer J, and symbols n and beta. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]] `):
print(`where  the actual size is denoted by a(n), and Av,2ndMom, ..., JthMom`):
print(`are the expressions, in terms of n and a(n), a(n+1), etc.`):
print(`for the average, 2ndMom (NOT about the mean), ..., all the way to the J-th moment`):
print(` (also not about the mean), followed by the generating function of a(n) `):
print(` Try: `):
print(` AvMomsAstraightE(1/(1-t-t^2*w),t,w,4,n,a); `):

elif nops([args])=1 and op(1,[args])=CheckAvMomsAM then
print(`CheckAvMomsAM(F,t,w,n0,J):  Checks AvMomsAM(F,t,w,n0,J): `):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs positive integers n0 and J. It returns a pair [L0,[Av,2ndMom, ..., JthMom]]]  with the sequence of numbers`):
print(` Try: `):
print(` CheckAvMomsAM(1/(1-t*(1+w)),t,w,20,4);`):

elif nops([args])=1 and op(1,[args])=FindOpe then
print(`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,`):
print(`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`):
print(`L1[n]=ope(n,N^(-1))L1[n] for n>=ord1+st1, for some st1<=MaxSt. It also returns the smallest st1. Try:`):
print(`FindOpe([seq(2^i,i=1..40)],[seq(i*2^i-2^(i-1),i=1..40)],n,N,0,1,5);`):

elif nops([args])=1 and op(1,[args])=FindOpe1 then
print(`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,`):
print(`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`):
print(`L1[n]=ope(n,N^(-1))L1[n] for n>=ord1+st1. Try:`):
print(`FindOpe1([seq(2^i,i=1..10)],[seq(i*2^i-2^(i-1),i=1..10)],n,N,0,1,1);`):

elif nops([args])=1 and op(1,[args])=LogCurveN then
print(`LogCurveN(F,t,w,N): inputs a rational function, F, of, t and w, let a(n,k) be the coeff. of w^k*t^n in the Maclaurin expansion of`):
print(`outputs the curve that is log(a(N,N*x))/N . Try: `):
print(`LogCurveN(1/(1-t*(1+w)),t,w,100);`):


elif nops([args])=1 and op(1,[args])=MDktw then
print(`MDktw(k,t,w): Inputs a positive integer k, 1<=k<=3, and outputs the pre-computed `):
print(`bi-variate rational function in t w, whose coefficient of t^n*w^r is the exact number`):
print(` of monomer-dimer tilings of a 2k by n strip with 2r monomers. Try: `):
print(`MDktw(2,t,w);`):
print(``):

elif nops([args])=1 and op(1,[args])=MomsSeq then
print(`MomsSeq(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,`):
print(`according to some statistics of the n-th entry in an infinite family of combinatorial objects indexed by n.`):
print(`It also inputs positive integers n0 and J. It returns a pair [L0,[L[1],...,L[J]]] `):
print(`where L0 is the list of the enumeration sequence from n=1 to n=n0, L1 is the list of numerators of the average, L2 is the list`):
print(`of numerators of the variance, etc. all the way to the J'th moment about the mean.`):
print(`Try: `):
print(`MomsSeq(1/(1-t*(1+w)),t,w,20,2);`):
print(`MomsSeq(MDktw(2,t,w),t,w,20,2);`):

elif nops([args])=1 and op(1,[args])=Pashet then
print(`Pashet(R,t,Pol): inputs a rational function R, in t, and a polynomial, Pol, in t, simplifies R modulo Pol. Try`):
print(`Pashet(t^3/(t^4+1),t,t^2-t-1);`):

elif nops([args])=1 and op(1,[args])=PW then
print(`PW(F,x,y,s,m): The Pemantle-Wilson approximate value for the coeff. of x^s*y^m in the rational function F(x,y). Try`):
print(`PW(1/(1-t-t^2*w,t,w,100,40); `):
print(`PW(1/(1-(1+w)*t),t,w,100,50); `):

elif nops([args])=1 and op(1,[args])=PWs then
print(`PWs(F,x,y,m,k): The Pemantle-Wilson done symbolically where s=m*k`):
print(`PWs(1/(1-t*(1+w)),t,w,m,k);`):

elif nops([args])=1 and op(1,[args])=PWsLL then
print(`PWsLL(F,x,y,k): The limit of log(PWs(F,x,y,m,k))/m as m goes to infinity. Try:`):
print(`PWsLL(1/(1-(1+w)*t),t,w,k); `):




elif nops([args])=1 and op(1,[args])=RatToAsy then
print(`RatToAsy(R,n,r): inputs a rational function in n and outputs its asymptotic expansion to order r. Try:`):
print(`RatToAsy((n+1)/(n+2),n,4); `):

elif nops([args])=1 and op(1,[args])=SeqT then
print(`SeqT(f,t,w,k,N): inputs a rational function f of t and w outputs the sequence of "diagonals with slope k"`):
print(`i.e. the sequence A(n*a,n*b), where A(n,m) is the coefficient of t^n*w^k in the Maclaurin expansion of f. `):
print(`Try: `):
print(`SeqT(1/(1-t-w*t^2),t,w,1/2,40); `):

elif nops([args])=1 and op(1,[args])=ShoreshG then
print(`ShoreshG(P,t): inputs a polynomial P in t, finds the largest root in floating-point. Try:`):
print(`Shoresh(t^3-t-1,t);`):

elif nops([args])=1 and op(1,[args])=Story then
print(`Story(F,t,w,J,n,beta): inputs a rational function F of t and w such that the coeff. of t^n, is the generating function`):
print(`according to some statistics defined on the n-th member of the family. It also inputs a positive integer J`):
print(`and symbols n and beta. It outputs an article about the average, variance, and the limits of the scaled higher moments `):
print(`(with exponetially small errors) up to the J's.`):
print(`For example for all ways of walking n stairs where at each step you either advance one or two steps, according`):
print(`to "number of 2s"`):
print(`Try: `):
print(`Story(1/(1-t-t^2*w),t,w,4,n,b);`):

elif nops([args])=1 and op(1,[args])=TestLLL then
print(`TestLLL(F,t,w,x,N1,N2): tests the local limit law for the coefficients of t^n*w^k in F with  k=av(n)+sd(n)*x for n from N1 to N2`):
print(`using asymptotic normality (assuming local limit rule). `):
print(`It returns a pair. The first is the appx. using the normal distibution`):
print(`and the second is how many standard-dviations it is from the mean. If this is very large or small the approximation is not good.`):
print(`Try: `):
print(`TestLLL(1/(1-t-t^2*w),t,w,0, 50,100);`):

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


#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:



#MDktw(k,t,w): Inputs a positive integer k, 1<=k<=3, and outputs the pre-computed 
#bi-variate rational function in t w, whose coefficient of t^n*w^r is the exact number
#of monomer-dimer tilings of a 2k by n strip with 2r monomers. Try:
#MDktw(2,t,w);
MDktw:=proc(k,t,w) local lu:

lu:=[-(t-1)/(-t^2*w+t^3-t*w-2*t+1), -(-3*t^5*w^2-2*t^6*w+3*t^4*w^2+t^7-3*t^5*w+3*t^3*w^2+5*t^4*w-3*t^2*w^2-5*t^5+11*t^3*w+2*t^4-9*t^2*w+6*t^3-2*t*w-3*t^2-2*t+1)/(-1+3*t
+5*t^2-15*t^3-9*t^4+5*t^6-9*t^7+21*t^5+t^9-t^8+5*t*w+13*t^2*w^2+t*w^2-27*t^3*w+14*t^5*w^3-3*t^4*w^4+24*t^5*w-41*t^4*w^2+4*t^3*w^3+3*t^5*w^4+6*t^6*w^2+29*t^5*w^2+t^3
*w^4-4*t^3*w^2-t^8*w^2-t^6*w^4-3*t^7*w^2-2*t^7*w^3-9*t^7*w+21*t^2*w+19*t^6*w-40*t^4*w+t^8*w-18*t^4*w^3+2*t^2*w^3), -(-1+8*t+8*t^2-208*t^3+234*t^4-4499*t^6-11612*t^7
+2172*t^5+31970*t^9+34832*t^8+424744*t^12-774090*t^14-153260*t^10+941714*t^16-759282*t^18+12*t*w+139*t^2*w^2+3*t*w^2-1103*t^3*w+65168*t^5*w^3-4295*t^4*w^4+24*t^3*w^
5+16575*t^5*w-1618*t^4*w^2-856*t^3*w^3+9*t^2*w^4+53778*t^5*w^4-55543*t^6*w^2+45980*t^5*w^2-115*t^3*w^4-1638*t^3*w^2+801157*t^8*w^2+1966598*t^8*w^4-346*t^6*w^8-\
283226*t^7*w^6+60*t^5*w^8-789870*t^7*w^4+7354*t^5*w^6-77490*t^6*w^4-396*t^4*w^6-391253*t^7*w^2+5*t^3*w^6+9103311*t^15*w^9+390098*t^20-108168*t^13+406*t^30-2525*t^28
+136*t^31-32*t^32-4*t^33+791368*t^19-112080*t^22+256782*t^23-551964*t^21+2735557*t^15*w-24866103*t^14*w^2-568326*t^7*w^5-20406*t^6*w^6+1064*t^5*w^7-86595130*t^18*w^
6+58473361*t^19*w^5+t^34+6972*t^24-79330*t^25+6776*t^26+15896*t^27-1968*t^29+37365655*t^19*w^3-79654653*t^18*w^4+1155*t^14*w^12-8908*t^13*w^13-280*t^12*w^14-5526330
*t^17*w^9+6659353*t^16*w^10+1385336*t^15*w^11+3161*t^10*w^12-6468775*t^13*w^9+451946*t^12*w^10+110344*t^11*w^11-26454*t^11-741912*t^17+5284981*t^19*w-21439711*t^18*
w^2+9280*t^16*w^14+80*t^13*w^15-3073060*t^14*w^10-724278*t^13*w^11-36029*t^12*w^12-5095010*t^17*w+41890001*t^16*w^8-68189138*t^17*w^5+423272*t^15-39246766*t^17*w^3+
112602589*t^16*w^4+40779621*t^15*w^5+1840373*t^11*w^9-41745*t^10*w^10-6904*t^9*w^11-93277559*t^14*w^6-14770698*t^13*w^7+11676092*t^12*w^8+359088*t^9*w-718980*t^7*w^
3-304*t^11*w^13+114774026*t^16*w^6+29833990*t^15*w^7-29401425*t^14*w^8+497349*t^16*w^12+78444*t^15*w^13+4060*t^14*w^14-121420*t^17*w^11+480*t^15*w^15+26296*t^5*w^5-\
108622*t^7*w-35060029*t^18*w^8+20041995*t^15*w^3-96653962*t^14*w^4-6245919*t^13*w^5+29073171*t^16*w^2-508932*t^11*w-4437744*t^10*w^2+3227448*t^9*w^3-99000*t^7*w^7-\
14117465*t^10*w^4+3260910*t^9*w^5+1331814*t^8*w^6-359670*t^13*w+13451971*t^12*w^2-6205243*t^11*w^3-3304*t^7*w^9-11790238*t^10*w^6+66786*t^9*w^7+237989*t^8*w^8+
451003*t^13*w^3+49309544*t^12*w^4-5888263*t^11*w^5-168080*t^9*w^9+7196*t^8*w^10+44995177*t^12*w^6+3296462*t^11*w^7-2472159*t^10*w^8-5*t^31*w^6-60*t^29*w^8+210*t^27*
w^10+420*t^25*w^12+75*t^23*w^14-9*t^32*w^4+16*t^30*w^6-214*t^28*w^8-3556*t^26*w^10-841*t^24*w^12+380*t^22*w^14+31*t^31*w^4-1090*t^29*w^6+7806*t^27*w^8+7285*t^25*w^
10-8869*t^23*w^12+280*t^21*w^14-31*t^32*w^2+355*t^30*w^4+674*t^28*w^6-36013*t^26*w^8+61509*t^24*w^10+2825*t^22*w^12-2120*t^20*w^14+286*t^31*w^2-7046*t^29*w^4+51170*
t^27*w^6-101310*t^25*w^8-122897*t^23*w^10+70778*t^21*w^12-6265*t^19*w^14+858*t^30*w^2-2310*t^28*w^4-98162*t^26*w^6+657875*t^24*w^8-521014*t^22*w^10+64769*t^20*w^12-\
6200*t^18*w^14-8500*t^29*w^2+150382*t^27*w^4-771052*t^25*w^6+623746*t^23*w^8+556371*t^21*w^10-177421*t^19*w^12+1400*t^17*w^14-9861*t^28*w^2-80770*t^26*w^4+1535046*t
^24*w^6-4605912*t^22*w^8+2337876*t^20*w^10-330277*t^18*w^12+115973*t^27*w^2-1531300*t^25*w^4+5312120*t^23*w^6-2478206*t^21*w^8-750813*t^19*w^10+37076*t^17*w^12+9305
*t^15*w^14+37543*t^26*w^2+1264385*t^24*w^4-10765441*t^22*w^6+16991757*t^20*w^8-5717385*t^18*w^10-873096*t^25*w^2+8636031*t^23*w^4-20212525*t^21*w^6+7760534*t^19*w^8
-1315020*t^17*w^10+393825*t^15*w^12+400*t^13*w^14+182528*t^24*w^2-8735052*t^22*w^4+39851459*t^20*w^6+3897672*t^23*w^2-28536123*t^21*w^4+43901515*t^19*w^6-15730668*t
^17*w^8+3856429*t^15*w^10-122158*t^13*w^12-75*t^11*w^14-2315755*t^22*w^2+34179438*t^20*w^4-10591279*t^21*w^2+55515336*t^19*w^4-55508292*t^17*w^6+18196666*t^15*w^8-\
2617399*t^13*w^10+8881*t^11*w^12+9579383*t^20*w^2+17618709*t^19*w^2-61296184*t^17*w^4+39247013*t^15*w^6-11512254*t^13*w^8+594449*t^11*w^10-420*t^9*w^12-17422960*t^
17*w^2+32814096*t^15*w^4-12808971*t^13*w^6+3418202*t^11*w^8-47285*t^9*w^10+8844251*t^15*w^2-713657*t^13*w^4-395488*t^11*w^6-283582*t^9*w^8-210*t^7*w^10-224285*t^13*
w^2-8430207*t^11*w^4+1402928*t^9*w^6-23318*t^7*w^8-2497980*t^11*w^2+4136976*t^9*w^4+1487736*t^9*w^2-3*t^33*w^2+36*t^31*w^5-184*t^29*w^7+1904*t^27*w^9+2984*t^25*w^11
-436*t^23*w^13-80*t^21*w^15+180*t^31*w^3-3380*t^29*w^5+22116*t^27*w^7-7160*t^25*w^9-48384*t^23*w^11+9908*t^21*w^13-480*t^19*w^15+197*t^31*w-8716*t^29*w^3+102470*t^
27*w^5-352554*t^25*w^7-63063*t^23*w^9+267106*t^21*w^11-41224*t^19*w^13-5329*t^29*w+158060*t^27*w^3-1233550*t^25*w^5+2467870*t^23*w^7+249601*t^21*w^9-506640*t^19*w^
11+13720*t^17*w^13+60530*t^27*w-1392464*t^25*w^3+7928789*t^23*w^5-9644166*t^21*w^7+782387*t^19*w^9-377572*t^25*w+6915613*t^23*w^3-28432563*t^21*w^5+23183242*t^19*w^
7+1428552*t^23*w-20650337*t^21*w^3-33834528*t^17*w^7-3429222*t^21*w+89*t^2*w-25422*t^6*w+415*t^4*w+250203*t^8*w-4447*t^4*w^3+66*t^2*w^3-420*t^26*w^11-6397918*t^14*w
-75889*t^6*w^3-49301752*t^18*w^3-30*t^4*w^7+27034151*t^12*w^7+131424941*t^16*w^5+677938*t^8*w^7+55084242*t^12*w^5-5749764*t^18*w-4184*t^6*w^7+77756660*t^16*w^7-\
110626824*t^14*w^5+54250*t^8*w^9-240*t^24*w^13-15066505*t^10*w^5-59826504*t^14*w^3-6578068*t^10*w^7+69196836*t^16*w^3-59613115*t^14*w^7+3329109*t^12*w^9-270*t^20*w^
15-10*t^22*w^15-552661*t^10*w^9+7543524*t^16*w-1219984*t^10*w-1920*t^4*w^5+1509743*t^8*w^3+3515685*t^12*w-52110*t^6*w^5-9601487*t^10*w^3+1882340*t^8*w^5+31438358*t^
12*w^3+3138*t^22*w^13-420*t^18*w^15+6000*t^24*w^11-1586*t^20*w^13+420*t^16*w^15+30*t^30*w^7-14266*t^26*w^9-75479*t^22*w^11-52524*t^18*w^13+270*t^14*w^15-248*t^28*w^
7+256319*t^24*w^9+519479*t^20*w^11+87484*t^16*w^13+10*t^12*w^15+200*t^30*w^5-66266*t^26*w^7-1892095*t^22*w^9-1578432*t^18*w^11+20506*t^14*w^13+6*t^32*w^3+2174*t^28*
w^5+1164960*t^24*w^7+7335535*t^20*w^9+2031660*t^16*w^11-5818*t^12*w^13+589*t^30*w^3-109012*t^26*w^5-8137433*t^22*w^7-15938821*t^18*w^9-499829*t^14*w^11+240*t^10*w^
13-43*t^32*w-7357*t^28*w^3+1598531*t^24*w^5+29686577*t^20*w^7+18316491*t^16*w^9-56671*t^12*w^11+1067*t^30*w-8405*t^26*w^3-10901014*t^22*w^5-61750572*t^18*w^7-\
11152997*t^14*w^9+11416*t^10*w^11-9710*t^28*w+681329*t^24*w^3+41756404*t^20*w^5+420*t^8*w^11+34511*t^26*w-5411154*t^22*w^3-94984535*t^18*w^5+29436*t^24*w+21408308*t
^20*w^3-694331*t^22*w+2722210*t^20*w)/(1-9*t-12*t^2+296*t^3-286*t^4+7268*t^6+26761*t^7-3896*t^5-103332*t^9-68735*t^8-1225424*t^12+2737322*t^14+365084*t^10-4163732*t
^16+4330900*t^18-18*t*w-241*t^2*w^2-8*t*w^2+1658*t^3*w-136527*t^5*w^3+13665*t^4*w^4-119*t^3*w^5-31058*t^5*w+6137*t^4*w^2+1485*t^3*w^3-37*t^2*w^4-t*w^3-118149*t^5*w^
4+45671*t^6*w^2-91026*t^5*w^2+155*t^3*w^4+2645*t^3*w^2-1490062*t^8*w^2-4100224*t^8*w^4-30802*t^6*w^8+2043729*t^7*w^6+670*t^5*w^8+2924085*t^7*w^4-14332*t^5*w^6-93623
*t^6*w^4+1762*t^4*w^6+1023966*t^7*w^2-37*t^3*w^6-150480235*t^15*w^9-3033924*t^20-70564*t^13-6048*t^30+20021*t^28-3300*t^31+794*t^32+168*t^33-3143968*t^19+1362430*t^
22-1771516*t^23+2898554*t^21-5077334*t^15*w+107548344*t^14*w^2+2825642*t^7*w^5-148721*t^6*w^6-305*t^5*w^7+786827772*t^18*w^6-598680250*t^19*w^5-48*t^34-335436*t^24+
723680*t^25+9824*t^26-194412*t^27+33023*t^29-264259094*t^19*w^3+596838548*t^18*w^4-329487*t^14*w^12+503706*t^13*w^13+34546*t^12*w^14+177311542*t^17*w^9-90721561*t^
16*w^10-32571185*t^15*w^11-42157*t^10*w^12+2940*t^9*w^13+70664285*t^13*w^9-9453064*t^12*w^10-2276228*t^11*w^11+208464*t^11+2145708*t^17-26641826*t^19*w+139600414*t^
18*w^2+525*t^11*w^15+16096*t^16*w^14-5013*t^13*w^15+38924791*t^14*w^10+12979743*t^13*w^11+542375*t^12*w^12+18320150*t^17*w-435866700*t^16*w^8+560213202*t^17*w^5-\
745654*t^15+211283862*t^17*w^3-720339076*t^16*w^4-248502359*t^15*w^5-15419662*t^11*w^9+1755287*t^10*w^10+164351*t^9*w^11+664542792*t^14*w^6+100523238*t^13*w^7-\
108706281*t^12*w^8-1183400*t^9*w+2158288*t^7*w^3-40937*t^11*w^13-1740*t^10*w^14-920777389*t^16*w^6-312953567*t^15*w^7+280182834*t^14*w^8-7970472*t^16*w^12-2829461*t
^15*w^13-227161*t^14*w^14+32893401*t^17*w^11-45809*t^15*w^15-2908*t^14*w^16-48*t^13*w^17-58530*t^5*w^5-3*t^3*w^7+262832*t^7*w+427193732*t^18*w^8-62832745*t^15*w^3+
524565059*t^14*w^4-10525496*t^13*w^5-3024*t^8*w^12-155083962*t^16*w^2+2819746*t^11*w+11648751*t^10*w^2-13420127*t^9*w^3+1101866*t^7*w^7+140*t^5*w^9+48785883*t^10*w^
4-22089108*t^9*w^5-4911343*t^8*w^6-2371306*t^13*w-46149151*t^12*w^2+37556808*t^11*w^3+114436*t^7*w^9-750*t^6*w^10+62142013*t^10*w^6-6078320*t^9*w^7-1996379*t^8*w^8-\
34354743*t^13*w^3-221254609*t^12*w^4+58557244*t^11*w^5+811360*t^9*w^9-193548*t^8*w^10+1950*t^7*w^11-281508133*t^12*w^6-4017086*t^11*w^7+23099098*t^10*w^8+9*t^31*w^
10+84*t^29*w^12+126*t^27*w^14+36*t^25*w^16+t^23*w^18+18*t^32*w^8-150*t^30*w^10-1092*t^28*w^12-564*t^26*w^14+114*t^24*w^16+10*t^22*w^18-158*t^31*w^8+2741*t^29*w^10+
4350*t^27*w^12-2819*t^25*w^14-456*t^23*w^16+45*t^21*w^18+210*t^32*w^6-266*t^30*w^8-16720*t^28*w^10+18845*t^26*w^12+11062*t^24*w^14-2692*t^22*w^16+120*t^20*w^18+121*
t^33*w^4-3056*t^31*w^6+22490*t^29*w^8-23059*t^27*w^10-117461*t^25*w^12+45965*t^23*w^14-3020*t^21*w^16+210*t^19*w^18+397*t^32*w^4+3023*t^30*w^6-84399*t^28*w^8+474909
*t^26*w^10-226391*t^24*w^12-39453*t^22*w^14+8636*t^20*w^16+252*t^18*w^18-10913*t^31*w^4+132583*t^29*w^6-541166*t^27*w^8-143967*t^25*w^10+1054265*t^23*w^12-345440*t^
21*w^14+33796*t^19*w^16+210*t^17*w^18+1661*t^32*w^2-5765*t^30*w^4-229223*t^28*w^6+2352346*t^26*w^8-5264634*t^24*w^10+2260185*t^22*w^12-364676*t^20*w^14+52208*t^18*w
^16+120*t^16*w^18-14314*t^31*w^2+335671*t^29*w^4-2553967*t^27*w^6+6001405*t^25*w^8+558191*t^23*w^10-3087228*t^21*w^12+492510*t^19*w^14+43948*t^17*w^16+45*t^15*w^18-\
30257*t^30*w^2+22430*t^28*w^4+3697857*t^26*w^6-22375271*t^24*w^8+28775761*t^22*w^10-10307250*t^20*w^12+1445010*t^18*w^14+18626*t^16*w^16+10*t^14*w^18+229506*t^29*w^
2-4232750*t^27*w^4+23878737*t^25*w^6-37085051*t^23*w^8+8751206*t^21*w^10-1349348*t^19*w^12+1148292*t^17*w^14+640*t^15*w^16+t^13*w^18+234636*t^28*w^2+1012491*t^26*w^
4-29422577*t^24*w^6+101644578*t^22*w^8-77911627*t^20*w^10+16227740*t^18*w^12-2063618*t^27*w^2+29039527*t^25*w^4-118255803*t^23*w^6+127466721*t^21*w^8-44135557*t^19*
w^10+11843130*t^17*w^12-493592*t^15*w^14-996*t^13*w^16-617173*t^26*w^2-14685395*t^24*w^4+141650850*t^22*w^6-264854884*t^20*w^8+112486410*t^18*w^10+11277982*t^25*w^2
-119384004*t^23*w^4+339567623*t^21*w^6-259018461*t^19*w^8+82189107*t^17*w^10-11046052*t^15*w^12+27345*t^13*w^14+36*t^11*w^16-2437087*t^24*w^2+89658539*t^22*w^4-\
421714181*t^20*w^6-38409024*t^23*w^2+302169976*t^21*w^4-582550428*t^19*w^6+322276735*t^17*w^8-77083682*t^15*w^10+3269433*t^13*w^12+607*t^11*w^14+23146042*t^22*w^2-\
300432012*t^20*w^4+81755204*t^21*w^2-464214848*t^19*w^4+589223474*t^17*w^6-241656439*t^15*w^8+35608061*t^13*w^10-431535*t^11*w^12+126*t^9*w^14-76425268*t^20*w^2-\
106389641*t^19*w^2+402650580*t^17*w^4-318718819*t^15*w^6+101758963*t^13*w^8-7415459*t^11*w^10+29550*t^9*w^12+77650285*t^17*w^2-144915366*t^15*w^4+53603671*t^13*w^6-\
18497443*t^11*w^8+527807*t^9*w^10+84*t^7*w^12-20129702*t^15*w^2-43453684*t^13*w^4+28956789*t^11*w^6-598728*t^9*w^8+19697*t^7*w^10-13443994*t^13*w^2+60351097*t^11*w^
4-15145013*t^9*w^6+428230*t^7*w^8+9*t^5*w^10+13954630*t^11*w^2-21125050*t^9*w^4-5351700*t^9*w^2-3*t^35+t^36+9*t^33*w^6+t^34*w^4+2*t^35*w^2-23*t^34*w^2+375*t^33*w^2-\
3*t^33*w^7-16*t^31*w^9+606*t^29*w^11+1092*t^27*w^13+45*t^25*w^15-60*t^23*w^17-t^35*w^3+9*t^33*w^5-993*t^31*w^7+8512*t^29*w^9+6007*t^27*w^11-26513*t^25*w^13+2131*t^
23*w^15+288*t^21*w^17-6*t^35*w+261*t^33*w^3-6596*t^31*w^5+57174*t^29*w^7-161004*t^27*w^9-275108*t^25*w^11+296436*t^23*w^13-52069*t^21*w^15+3960*t^19*w^17+284*t^33*w
-15761*t^31*w^3+244070*t^29*w^5-1324760*t^27*w^7+1395868*t^25*w^9+2060625*t^23*w^11-1333629*t^21*w^13+166665*t^19*w^15+4464*t^17*w^17-7776*t^31*w+326016*t^29*w^3-\
3813048*t^27*w^5+14290038*t^25*w^7-10152931*t^23*w^9-2826897*t^21*w^11+669466*t^19*w^13+266837*t^17*w^15+612*t^15*w^17+111392*t^29*w-3458309*t^27*w^3+30004484*t^25*
w^5-78330452*t^23*w^7+47692357*t^21*w^9-11708193*t^19*w^11+3929700*t^17*w^13-883350*t^27*w+21217854*t^25*w^3-135585322*t^23*w^5+239658681*t^21*w^7-120862268*t^19*w^
9+4160854*t^25*w-79543367*t^23*w^3+367406801*t^21*w^5-438362897*t^19*w^7-12193656*t^23*w+185485301*t^21*w^3+484882525*t^17*w^7+22711950*t^21*w-131*t^2*w+34393*t^6*w
+280*t^4*w-477162*t^8*w+13779*t^4*w^3-152*t^2*w^3+126726*t^26*w^11+24946198*t^14*w-7172*t^6*w^3+342775504*t^18*w^3+191*t^4*w^7-203586598*t^12*w^7-936771324*t^16*w^5
-3687284*t^8*w^7-290186428*t^12*w^5+35683278*t^18*w-86435*t^6*w^7-707611257*t^16*w^7+685133706*t^14*w^5-755540*t^8*w^9+6604*t^24*w^13+63470229*t^10*w^5+287879712*t^
14*w^3+44781755*t^10*w^7-402210356*t^16*w^3+491014073*t^14*w^7-40998909*t^12*w^9-8975*t^20*w^15-18579*t^22*w^15+8150001*t^10*w^9-36958030*t^16*w+3003282*t^10*w+6895
*t^4*w^5-2835964*t^8*w^3-10925802*t^12*w-156675*t^6*w^5+28238197*t^10*w^3-4966381*t^8*w^5-121971339*t^12*w^3+231940*t^22*w^13+323228*t^18*w^15-1476209*t^24*w^11-\
2468634*t^20*w^13+70605*t^16*w^15+1147*t^30*w^7+1225703*t^26*w^9+9890773*t^22*w^11+5189844*t^18*w^13-37537*t^14*w^15-159136*t^28*w^7-12752169*t^24*w^9-31760245*t^20
*w^11-1275484*t^16*w^13+2641*t^12*w^15+1107*t^30*w^5+3432661*t^26*w^7+61568761*t^22*w^9+45508458*t^18*w^11-700318*t^14*w^13+579*t^32*w^3-185885*t^28*w^5-29284850*t^
24*w^7-157581774*t^20*w^9-30708649*t^16*w^11+203142*t^12*w^13-20414*t^30*w^3+2694571*t^26*w^5+133371575*t^22*w^7+238516436*t^18*w^9+7205905*t^14*w^11-14005*t^10*w^
13+1698*t^32*w+206888*t^28*w^3-23383864*t^24*w^5-368197833*t^20*w^7-218994974*t^16*w^9-429297*t^12*w^11-23711*t^30*w-303265*t^26*w^3+123635080*t^22*w^5+639455064*t^
18*w^7+122591935*t^14*w^9+128802*t^10*w^11+143410*t^28*w-6986027*t^24*w^3-396207778*t^20*w^5-31832*t^8*w^11-289646*t^26*w+52430694*t^22*w^3+776882916*t^18*w^5-\
1010546*t^24*w-177219840*t^20*w^3+7722658*t^22*w-21849264*t^20*w-t^32*w^9-t^4*w^9-9*t^24*w^17-84*t^26*w^15-102*t^22*w^17-126*t^28*w^13+2097*t^24*w^15+1737*t^20*w^17
-36*t^30*w^11-37*t^26*w^13+5292*t^18*w^17-5288*t^28*w^11+2313*t^16*w^17-426*t^30*w^9-6*t^14*w^17+47*t^32*w^7-40136*t^28*w^9-9*t^12*w^17-3*t^34*w^5+333*t^32*w^5-24*t
^34*w^3-84*t^10*w^15-41*t^34*w-126*t^8*w^13-36*t^6*w^11-6522*t^6*w^9-3*t^2*w^5)]:
lu[k]:
end:


#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:


#AsySizeT(R,t,beta): inputs a rational function R of the variable t, that is the  generating function  of some sequence of positive
#integers a(n). Returns [A,beta,PolSatisfyingbeta] where the 
# 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). The output is the triple [A,beta,POLinBeta].  Try:
#AsySizeT(1/(1-t-t^2),t,beta);
AsySizeT:=proc(R,t,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):
bot1:=subs(t=beta,bot1):

[A, beta,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:



#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:




#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:

#AsyData(F,t,w,beta): inputs a rational function F in the variables t and w such that the coeff. of t^n is the generating function
#according to some statistic of the n-th member of an sequence of combinatorial objects, outputs the triple
#[A,c1,c2,beta,pol] where c1 and c2 are expressions in terms of the algebraric number beta, (the largest positive root
#of the polynomial in beta, and the meaning of c1 and c2 are that the asymptotics of the average of the statistic is c1*n and 
#the asymptotics of the variance is c2*n. The asymptotics for the cardinality of the n-th family is A*beta^n.  Try:
#AsyData(1/(1-t-w*t^2),t,w,beta);
AsyData:=proc(F,t,w,beta) local c1,c2,n,ku,A,bot1:

 ku:=AsySizeT(subs(w=1,F),t,beta):
A:=ku[1]:
bot1:=ku[3]:

 ku:=AvMomsA(F,t,w,2,n,beta)[2]:

c1:=coeff(ku[1],n,1):
c2:=coeff(ku[2],n,1):

[A,c1,c2,beta,bot1]:

end:

#AsyDataFl(F,t,w): Like AsyDataFl(F,t,w,beta) (q.v.) but in floating-point. 
#It returns the four numbers [A,beta,c1,c2] such that the cardinality of the n-th member in the sequence is 
#asympototically, A*beta^n, the expectation of the statistic is asymptotically c1*n, and the variance is asymptotically c2*n.
#Try:
#AsyDataFl(1/(1-t-w*t^2),t,w);
AsyDataFl:=proc(F,t,w) local gu,beta,bot1,b1,A,c1,c2:

gu:=AsyData(F,t,w,beta):

bot1:=gu[5]:

b1:=ShoreshG(bot1,beta):

if abs(coeff(b1,I,1))>10^(-5) then
 RETURN(FAIL):
else
 b1:=coeff(b1,I,0):
fi:

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

if abs(coeff(A,I,1))>10^(-5) then
 RETURN(FAIL):
else
 A:=coeff(A,I,0):
fi:

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

if abs(coeff(c1,I,1))>10^(-5) then
 RETURN(FAIL):
else
 c1:=coeff(c1,I,0):
fi:

c2:=subs(beta=b1,gu[3]):

if abs(coeff(c2,I,1))>10^(-5) then
 RETURN(FAIL):
else
 c2:=coeff(c2,I,0):
fi:

[A,b1,c1,c2]:


end:



#Pashet(R,t,Pol): inputs a rational function R, in t, and a polynomial, Pol, in t, simplifies R modulo Pol. Try
#Pashet(t^3/(t^4+1),t,t^2-t-1);
Pashet:=proc(R,t,Pol) :
 normal(rem(numer(R),Pol,t)/rem(denom(R),Pol,t)):
end:




#Alpha(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,Variance,3rdStandarizedMoment, ..., JthStandarized]]] 
#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, Variance, Skewness^2, Kurtosis..., JthAlphaCoeffi. are the asymptotic expressions
#in terms of beta  and n for the average, variance ..., all the way to the stadarized J-th moment, BUT 
#Note that the odd-standarized moments are squared.
#(also about the mean), followed by the minimal polynomial in t, such that beta is its largest root.
#Try:
#Alpha(1/(1-t-t^2*w),t,w,4,n,beta);
Alpha:=proc(F,t,w,J,n,beta) local  gu,mu,m2,lu,ku,pol,j:

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

pol:=gu[3]:


mu:=gu[2]:


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


m2:=mu[2]:

for j from 3 to J do

if j mod 2=0 then
 ku:=normal(mu[j]/m2^(j/2)):
else
 ku:=normal(mu[j]^2/m2^j):
fi:
  ku:=Pashet(ku,beta,pol):
 lu:=[op(lu),ku]:
od:

[gu[1],lu,gu[3],gu[4]]:

end:




#AlphaFl(F,t,w,J,n): Like Alpha(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 STANDARIZED moment
#Try:
#AlphaFl(1/(1-t-t^2*w),t,w,4,n);
AlphaFl:=proc(F,t,w,J,n) local beta,gu,b1:

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

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



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

#RatToAsy(R,n,r): inputs a rational function in n and outputs its asymptotic expansion to order r. Try:
#RatToAsy((n+1)/(n+2),n,r):
RatToAsy:=proc(R,n,r) local gu,i,x:

if degree(numer(R),n)>degree(denom(R),n) then
 print(`Top degree must be <=0  bottom degree `):
 RETURN(FAIL):
fi:

gu:=normal(subs(n=1/x,R)):
gu:=taylor(gu,x=0,r+1):
add(coeff(gu,x,i)/n^i,i=0..r):
end:




#AlphaAsy(F,t,w,J,n,beta,r): Like AlphaAsy(F,t,w,J,n,beta) (q.v.) but with an asymptotic expression in n to order r
#of the standarized 3rd- to J-th moments (for the odd ones its square).
#Try:
#AlphaAsy(1/(1-t-t^2*w),t,w,6,n,beta,2);
AlphaAsy:=proc(F,t,w,J,n,beta,r) local  gu,j,ku,lu,mu,pol,i:

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

pol:=gu[3]:


mu:=gu[2]:


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

for j from 3 to J do

 ku:=mu[j]:

 ku:=RatToAsy(ku,n,r):
 ku:=add( Pashet(coeff(ku,n,-i),beta, pol)*n^(-i),i=0..r):
 lu:=[op(lu),ku]:
od:

[gu[1],lu,gu[3],gu[4]]:
end:


#AppxAnk(F,t,w,n1,k1): the approximation of the coeff. of t^n1*w^k1 in the Taylor expansion of the rational function F in t,w
#using asymptotic normality (assuming local limit rule). 
#It returns a pair. The first is the appx. using the normal distibution
#and the second is how many standard-dviations it is from the mean. If this is very large or small the approximation is not good.
#Try:
#AppxAnk(1/(1-t-t^2*w),t,w,100,50);
AppxAnk:=proc(F,t,w,n1,k1) local gu,val,av1,sd1,n,kama:
gu:=AvMomsAfl(F,t,w,2,n):
val:=subs(n=n1,gu[1]):
av1:=subs(n=n1,gu[2][1]):
sd1:=sqrt(subs(n=n1,gu[2][2])):

kama:=(k1-av1)/sd1:
[evalf(val/sqrt(2*Pi)/sd1*exp(-kama^2/2)),kama]:

end:


#ExactAnk(F,t,w,n1,k1): the Exact value of the coeff. of t^n1*w^k1 in the Taylor expansion of the rational function F in t,w
#using asymptotic normality (assuming local limit rule). Try:
#ExactAnk(1/(1-t-t^2*w),t,w,100,50);
ExactAnk:=proc(F,t,w,n1,k1):
coeff(expand(coeff(taylor(F,t=0,n1+1),t,n1)),w,k1):
end:


#TestLLL(F,t,w,x,N1,N2): tests the local limit law for the coefficients of t^n*w^k in F with  k=av(n)+sd(n)*x for n from N1 to N2
#using asymptotic normality (assuming local limit rule). 
#It returns a pair. The first is the appx. using the normal distibution
#and the second is how many standard-dviations it is from the mean. If this is very large or small the approximation is not good.
#Try:
#TestLLL(1/(1-t-t^2*w),t,w,50,100);
TestLLL:=proc(F,t,w,x,N1,N2) local gu,val,av1,sd1,n,kama,k1,lu,n1:
gu:=AvMomsAfl(F,t,w,2,n):

lu:=[]:

for n1 from N1 to N2 do
val:=subs(n=n1,gu[1]):
av1:=subs(n=n1,gu[2][1]):
sd1:=sqrt(subs(n=n1,gu[2][2])):
k1:=trunc(av1+sd1*x):
lu:=[op(lu),evalf(val/sqrt(2*Pi)/sd1*exp(-x^2/2)/ExactAnk(F,t,w,n1,k1))]:
od:

gu: 

end:





#LogCurveN(F,t,w,N): inputs a rational function, F, of, t and w, let a(n,k) be the coeff. of w^k*t^n in the Maclaurin expansion of
#outputs the curve that is log(a(N,N*x))/N . Try 
#LogCurveN(1/(1-t*(1+w)),t,w,100);
LogCurveN:=proc(F,t,w,N) local gu,i:
gu:=expand(coeff(taylor(F,t=0,N+1),t,N)):

[seq([i/N,evalf(log(coeff(gu,w,i))/N)],i=0..degree(gu,w))]:

end:



#QH:=proc(H,x,y):
#-x*diff(H,x)*(y*diff(H,y))^2-y*diff(H,y)*(x*diff(H,x))^2-
#((y*diff(H,y))^2*(x^2*diff(H,x,x))+ (x*diff(H,x))^2*(y^2*diff(H,y,y))
#-2*x*diff(H,x)*y*diff(H,y)*x*y*diff(H,x,y)):
#end:

#PW(F,x,y,s,m): The Pemantle-Wilson approximate value for the coeff. of x^s*y^m in the rational function F(x,y). Try
#PW(1/(1-t-t^2*w,t,w,100,40),
PW:=proc(F,x,y,s,m) local F1,G,H,x0,y0,Q,lu:

F1:=normal(F):

G:=numer(F1): H:=denom(F1):

lu:=fsolve({H, m*x*diff(H,x)-s*y*diff(H,y)},{x,y}):
x0:=subs(lu,x):
y0:=subs(lu,y):

Q:=-x*diff(H,x)*(y*diff(H,y))^2-y*diff(H,y)*(x*diff(H,x))^2-
((y*diff(H,y))^2*(x^2*diff(H,x,x))+ (x*diff(H,x))^2*(y^2*diff(H,y,y))
-2*x*diff(H,x)*y*diff(H,y)*x*y*diff(H,x,y)):
Q:=normal(Q):

abs(evalf(1/sqrt(2*Pi)*subs({x=x0,y=y0},G)*x0^(-s)*y0^(-m)*sqrt(-y0*subs({x=x0,y=y0},diff(H,y))/m/subs({x=x0,y=y0},Q)))):
end:


#PWs(F,x,y,m,k): The Pemantle-Wilson done symbolically where s=m*k
#PWs(1/(1-t-t^2*w,t,w,m,k),
PWs:=proc(F,x,y,m,k) local F1,G,H,x0,y0,Q,lu:

F1:=normal(F):

G:=numer(F1): H:=denom(F1):

lu:=solve({H, x*diff(H,x)-k*y*diff(H,y)},{x,y}):

x0:=subs(lu,x):
y0:=subs(lu,y):

#RETURN([x0,y0]):

Q:=-x*diff(H,x)*(y*diff(H,y))^2-y*diff(H,y)*(x*diff(H,x))^2-
((y*diff(H,y))^2*(x^2*diff(H,x,x))+ (x*diff(H,x))^2*(y^2*diff(H,y,y))
-2*x*diff(H,x)*y*diff(H,y)*x*y*diff(H,x,y)):
Q:=normal(Q):

1/sqrt(2*Pi)*subs({x=x0,y=y0},G)*x0^(-m*k)*y0^(-m)*sqrt(-y0*subs({x=x0,y=y0},diff(H,y))/m/subs({x=x0,y=y0},Q));
normal(%):
end:


#PWsLL(F,x,y,k): The limit of log(PWs(F,x,y,m,k))/m as m goes to infinity. Try:
#PWsLL(1/(1-(1+w)*t),t,w,k),
PWsLL:=proc(F,x,y,k) local m, F1,G,H,x0,y0,lu:

F1:=normal(F):

G:=numer(F1): H:=denom(F1):

lu:=solve({H, x*diff(H,x)-k*y*diff(H,y)},{x,y}):

x0:=subs(lu,x):
y0:=subs(lu,y):


-k*log(x0)-log(y0):
end:



#Story(F,t,w,J,n,beta): inputs a rational function F of t and w such that the coeff. of t^n, is the generating function
#according to some statistics defined on the n-th member of the family. It also inputs a positive integer J
#and symbols n and beta. It outputs an article about the average, variance, and the limits of the scaled higher moments 
#(with expontially small errors) up to the J's.
#For example for all ways of walking n stairs where at each step you either advance one or two steps, according
#to "number of 2s"
#Try:
#Story(1/(1-t-t^2*w),t,w,4,n,b);
Story:=proc(F,t,w,J,n,beta) local gu,P,bf,j:
print(``):
print(`The Average, Variance, and Scaled Moments Up to the`, J, `-th of a Certain Combinatorial Random Variable`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Suppose that there is an infinite family of cominatorial families, on which there is a combinatorial random variable`):
print(` Let `, P[n](w) , `be the weight-enumerator according to that random variable of the n-th member of our family. Suppose that`):
print(`you found out (e.g. by the transfer matrix method) that the bi-variate generating function, `, Sum(P[n](w)*t^n,n=0..infinity), ` equals `, F):
print(``):
print(`In other words`):
print(``):
print(Sum(P[n](w)*t^n,n=0..infinity)= F):
print(``):

##
gu:=AvMomsAE(F,t,w,J,n,a):
print(`Definition: Let a(n) be defined by the ordinary generating function`):
print(Sum(a(n)*t^n,n=0..infinity)=gu[2]):
gu:=gu[1]:
print(``):
print(`then, of course, a(n), is the number of elements in the n-th family.`):
print(``):
print(`The EXACT expression, in terms of a(n) and n, for the EXPECTATION is`):
print(``):
print(normal(gu[1]/a(n))):
print(``):
print(`and in Maple notation`):
print(``):
lprint(normal(gu[1]/a(n))):
print(``):
if J>1 then
 print(`The EXACT expression, in terms of a(n) and n, for the VARIANCE is`):
print(``):
print(normal(gu[2]/a(n)^2)):
print(``):
print(`and in Maple notation`):
lprint(normal(gu[2]/a(n)^2)):
print(``):
fi:

for j from 3 to J do
 print(`The EXACT expression, in terms of a(n) and n, for the`, j, `-th moment about the mean is`):
print(``):
print(normal(gu[j]/a(n)^j)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(normal(gu[j]/a(n)^j)):
od:
print(``):


gu:=Alpha(F,t,w,J,n,beta):
print(``):
print(` Let `, beta, `be the largest positive root  of the polynomial equation`):
print(``):
print(gu[3]=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[3]=0):
print(``):
print(`whose floating-point approximation is`):
print(``):
bf:=ShoreshG(gu[3],beta):
print(bf):
print(``):
print(`Then the size of the n-th family (i.e. straight enumeration) is very close to `):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):
print(`and in floating point`):
print(``):
print(subs(beta=bf,gu[1])):
print(``):
print(`The average of the statistics is, asymptotically `):
print(``):
print(collect(gu[2][1],n)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(collect(gu[2][1],n)):
print(``):
print(`and in floating-point`):
print(``):
lprint(subs(b=bf,collect(gu[2][1],n))):
print(``):

print(`The variance of the statistics is, asymptotically `):
print(``):
print(collect(gu[2][2],n)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(collect(gu[2][2],n)):
print(``):
print(`and in floating-point`):
print(``):
lprint(subs(b=bf,collect(gu[2][2],n))):
print(``):

if J>=3 then

print(`The skewness of the statistics is, asymptotically `):
print(``):
print(collect(gu[2][3],n)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(collect(gu[2][3],n)):
print(``):
print(`and in floating-point`):
print(``):
lprint(subs(b=bf,collect(gu[2][3],n))):
print(``):
fi:
	

if J>=4 then

print(`The kurtosis of the statistics is, asymptotically `):
print(``):
print(collect(gu[2][4],n)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(collect(gu[2][4],n)):
print(``):
print(`and in floating-point`):
print(``):
lprint(subs(b=bf,collect(gu[2][4],n))):
print(``):

fi:
	
for j from 5 to J do
print(`The standardized`, j, `-th moment (about the mean) of the statistics is, asymptotically `):
print(``):
print(collect(gu[2][j],n)):
print(``):
print(`and in Maple notation`):
print(``):
lprint(collect(gu[2][j],n)):
print(``):
print(`and in floating-point`):
print(``):
lprint(subs(b=bf,collect(gu[2][j],n))):
print(``):
od:

gu:=AlphaAsy(F,t,w,6,n,beta,2):

print(`Finally here is the asymptotic expressions, to order 2,  of the standarized third to`, J, `-th moment`):
	       
for j from 3 to J do
print(` The `, j, `-th standardized moment is`, gu[2][j]):
print(``):
print(`and in floating-point`, subs(beta=bf, gu[2][j])):
print(``):
od:

end:


ez:=proc(): print(` SeqT(f,t,w,k,N), fmkN(f,t,w,k,N), fmkNdata(f,t,w,A,N),  Zinn(resh) `): end:

sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

if resh[-1]=resh[-2] then
  RETURN(FAIL):
fi:
 
if resh[nops(resh)]=0 or resh[nops(resh)-1]=0 then
 RETURN(FAIL):
fi:
n1:=nops(resh)-1:
s1:=evalf(sn(resh,n1)):
s2:=evalf(sn(resh,n1-1)):

if abs(s1)<10^(-5)  or abs(s2)<10^(-5) or abs(s1-s2)<10^(-5) then
 RETURN(FAIL):
fi:

theta:=evalf(2*(s1+s2)/(s1-s2)^2):

mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end:
 

#SeqT(f,t,w,k,N): inputs a rational function f of t and w outputs the sequence of "diagonals with slope k"
#i.e. the sequence A(n*a,n*b), where A(n,m) is the coefficient of t^n*w^k in the Maclaurin expansion of f. 
#Try:
#SeqT(1/(1-t-w*t^2),t,w,1/2,40);

SeqT:=proc(f,t,w,k,N) local f1,a,b,i,n,L:
option remember:
a:=numer(k):
b:=denom(k):

f1:=taylor(f,t=0,N*b+1):

L:=[seq(coeff(f1,t,i),i=1..N*b)]:



[seq(coeff(L[b*n],w,n*a),n=1..nops(L)/b)]:


end:

#fmkN(f,t,w,k,N): Let A(k,n) (k rational, k=a/b, say, with a and b integers), be the coefficient of t^(n*b)*w^(a*n)
#finds (numerically) the limit of log(A(k,n))/(n*b) as n goes to infinity
#Try:
#fmkN(1/(1-t-t*w^2),1/2,100);
fmkN:=proc(f,t,w,k,N) local L,gu,b:
b:=denom(k):

L:=SeqT(f,t,w,k,N):

gu:=Zinn(L):


if gu=FAIL then
 RETURN(FAIL):
else
 gu:=gu[2]:
fi:

log(gu)/b:
end:


#fmkNdata(f,t,w,A,N): a table of [k,fmkN(f,t,w,k,N)] for all k=a/b with gcd(a,b)=1, 1<=a<b<=A. N is a fairly large positive integer.
#The higher N, the better the appxroximation (but the longer it takes)
#Try:
#fmkNdata(1/(1-t-t^2*w),t,w,3,100);
fmkNdata:=proc(f,t,w,A,N) local a,b,gu,T,i,lu:

gu:=[]:

for b from 2 to A do
for a from 1 to b-1 do
if gcd(a,b)=1 then 
 lu:=fmkN(f,t,w,a/b,N):
 if lu<>FAIL then
T[a/b]:=lu:
 gu:=[op(gu),a/b]:
fi:
fi:
od:
od:

gu:=sort(gu):
[seq([gu[i],T[gu[i]]],i=1..nops(gu))]:

end:






#LogCurveN(F,t,w,N): inputs a rational function, F, of, t and w, let a(n,k) be the coeff. of w^k*t^n in the Maclaurin expansion of
#outputs the curve that is log(a(N,N*x))/N . Try 
#LogCurveN(1/(1-t*(1+w)),t,w,100);
LogCurveN:=proc(F,t,w,N) local gu,i:
gu:=expand(coeff(taylor(F,t=0,N+1),t,N)):

[seq([i/N,evalf(log(coeff(gu,w,i))/N)],i=0..degree(gu,w))]:

end:



#AvMomsAstraightE(F,t,w,J,n,a): 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  the actual size is denoted by a(n), and Av,2ndMom, ..., JthMom
#are the expressions, in terms of n and a(n), a(n+1), etc.
#for the average, 2ndMom (NOT about the mean), ..., all the way to the J-th moment
#(also not about the mean), followed by the generating function of a(n)
#Try:
#AvMomsAstraightE(1/(1-t-t^2*w),t,w,4,n,a);
AvMomsAstraightE:=proc(F,t,w,J,n,a) local ord1,bot1,mu,lu,n0,N,j,ope1,i1,ku,beta:
ku:=AsySize( normal(subs(w=1,F)),t,n,beta):
bot1:=ku[2]:

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

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


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:=add(coeff(ope1,N,-i1)*a(n-i1),i1=0..-ldegree(ope1,N)):

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

[lu,subs(w=1,F)]:

end:


#AvMomsAM(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]]]  with the sequence of numbers
#for the moments about the mean
#AvMomsAM(1/(1-t*(1+w)),t,w,10,2);
#AvMomsAM(MDktw(2,t,w),t,w,10,2);
AvMomsAM:=proc(F,t,w,n0,J) local gu,i,k,r,N0:
gu:=AvMoms(F,t,w,n0,J):
N0:=gu[1]:
gu:=gu[2]:

[N0,
[gu[1],seq([seq(add((-1)^r*binomial(k,r)*gu[k-r][i]*N0[i]^(k-r-1)*gu[1][i]^r,r=0..k-1)+(-1)^k*gu[1][i]^k,i=1..n0)],k=2..J)]]:
end:


#CheckAvMomsAM(F,t,w,n0,J):  Checks AvMomsAM(F,t,w,n0,J): 
#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]]]  with the sequence of numbers
#Try:
#CheckAvMomsAM(1/(1-t*(1+w)),t,w,20,4);

CheckAvMomsAM:=proc(F,t,w,n0,J) local mu,gu,i,mu0,j,r:
mu:=AvMomsAM(F,t,w,n0,J):

mu0:=mu[1]:
mu:=mu[2]:

mu:=[seq([seq(mu[r][i]/mu0[i]^r,i=1..n0)],r=1..J)]:
gu:=[seq(expand(coeff(taylor(F,t=0,n0+1),t,i)),i=1..n0)]:

gu:=[seq(AveAndMoms(gu[i],w,J),i=1..n0)]:

gu:=[seq([seq(gu[i][j],i=1..n0)],j=1..J)]:

evalb(gu=mu):

end:






#AvMomsAE(F,t,w,J,n,a): 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  the actual size is denoted by a(n), and Av,2ndMom, ..., JthMom
#are the expressions, in terms of n and a(n), a(n+1), etc.
#for the average, and the moments ABOUT THE MEAN, ..., all the way to the J-th moment
#(also about the mean), followed by the generating function of a(n)
#Try:
#AvMomsAE(1/(1-t-t^2*w),t,w,4,n,a);
AvMomsAE:=proc(F,t,w,J,n,a) local gu,mu1,mu,av,r,akha,k:
gu:=AvMomsAstraightE(F,t,w,J,n,a):
akha:=gu[2]:
gu:=gu[1]:
mu:=[gu[1]]:
av:=gu[1]:
for k from 2 to J do
 mu1:=expand((-av)^k+add((-av)^r*binomial(k,r)*gu[k-r]*a(n)^(k-r-1),r=0..k-1)):
 mu:=[op(mu),mu1]:
od:
[mu,akha]:
end: