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


print(`Written:  Nov. 2021-Jan. 2022: tested for Maple 2020 `):
with(combinat):


print():
print(`This is  ETSIM.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad  and Doron Zeilberger's  article: `):
print(` Automated Counting and Statistical Analysis of Labeled Trees with Degree Restrictions `):

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

print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print():
 
print(`-----------------------------`):

print(`For a list of the supporting procedures type ezra1()`):
print(`For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():

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

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

print(`For a list of the Statistics procedures type ezraStat()`):
print(`For specific help type "ezra(procedure_name);" `):
print():

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


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

print(`For a list of the pre-derrived formulas for the average, variance, covariance etc., type ezraF(); `):
print(`For specific help type "ezra(procedure_name);" `):
print():

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

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

print(`For a list of the Findrec procedures type ezraFindrec()`):
print(`For specific help type "ezraFindrec(procedure_name);" `):

print():

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





ezraF:=proc()
if args=NULL then

print(` The pre-derived procedures with formulas for the average, variance, covariance of the random variables, X_d:= "number of vertices of degree d" `):
print(`LimAveFor, LimCorFor, LimCovFor, LimVarFor `):

else
ezra(args):

fi:
end:


ezraC:=proc()
if args=NULL then


print(`The CHECKING  procedures are`):
print(` CheckLTseq, CheckLTseqF, CheckKTseqF, CheckLTseqWE,  CheckLTseqFwe, CheckLTmom1,  CheckLTmom1am, DegSeq, Wt`):
else
ezra(args):

fi:
end:

ezraStat:=proc()
if args=NULL then

print(`The Statistics  procedures are`):
print(` Fc, LimMixedMoms, LimNor, LTrecG, LTmom, LTmom1am, LTmom1amL, LTmom2am, LTmom2amL, LTseqWE , LTseqFwe, LTstat, LTstatJ1, MixMomNorF,  ProveAsyNor, ProveJointAsyNor, ScaledMixedMom`):

else
ezra(args):

fi:
end:

ezraS:=proc()
if args=NULL then

print(`The Story  procedures are`):
print(`  InfoLTv, InfoLTeven, LTFpaper, LTpaper, LTpaperPlain , StatPaper  `):

else
ezra(args):

fi:
end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  As, CG, FEtoSeq,  GG, GGsp, GP, Kids, Kids1, LTmom1, LTmom2, LTrec, LTseqEven, LTseqFwe1, LTseqSlow, RG,  STd, u, STfNu, STmtt, STmttF, STmttP, STnu, SubsFPS,  T1rn, T2rn,T3rn, Zinn `):

else
ezra(args):

fi:
end:

ezra:=proc()
if args=NULL then

print(`  SPAN.txt: A Maple package to generate, randomly-generate, and count, spanning trees with degree restrictions `):

print(`The MAIN procedures are`):
print(` InfoLT, LTseq, LTseqF, ST, STf, STp, `):


elif nargs=1 and args[1]=As then
print(`As(ope,n,N,k): Outputs the connective constant and criticacl exponent. Try`):
print(`As(-(n-1)*(n+2)*(n+1)/(n+4)-(n+2)*(2*n+3)/(n+4)*(n+1)*N+(n+1)*(n+2)*N^2,n,N)`):
print(`As(-1/2*(n-1)*(n+1)*(n+3)^2*n/(n+4)/(2*n+7)-1/2*(n+3)*n/(n+4)/(2*n+7)*(n+1)*N-1/4*(n+3)*(21*n^3+147*n^2+326*n+224)/(n+4)/(2*n+7)*(n+1)*N^2+(n+1)*(n+2)*(n+3)*N^3,n,N)`):


elif nargs=1 and args[1]=CG then
print(`CG(n): inputs positive integer n, and outputs the complete graph on {1,...n}. Try:`);
print(`CG(5);`):


elif nargs=1 and args[1]=CheckLTmom1 then
print(`CheckLTmom1(d,N,a): checks LTmom1(d,n,a) for n=1..N. Try:`):
print(`CheckLTmom1(1,6,1);`):


elif nargs=1 and args[1]=CheckLTmom1am then
print(`CheckLTmom1am(d,N,a): checks LTmom1am(d,n,a) for n=1..N. Try:`):
print(`CheckLTmom1am(1,6,2);`):

elif nargs=1 and args[1]=CheckLTseq then
print(`CheckLTseq: Checks empirically LTseq. Try:`):
print(` CheckLTseq({1,2,3},6);`):


elif nargs=1 and args[1]=CheckLTseqF then
print(`CheckLTseqF: Checks empirically LTseqF. Try:`):
print(` CheckLTseqF({2,3},6);`):

elif nargs=1 and args[1]=CheckLTseqWE then
print(`CheckLTseqWE(P,N): Checks LTseqWE(P,N,g). Try:`):
print(` CheckLTseqWE([1,2,3],5); `):


elif nargs=1 and args[1]=CheckLTseqFwe then
print(`CheckLTseqFwe(P,N): Checks LTseqWE(P,N,g). Try:`):
print(`CheckLTseqFwe([1,2,3],6);`):

elif nargs=1 and args[1]=DegSeq then
print(`DegSeq(G,n) : Given a graph G on n vertices, outputs the list of length n such that its i-th entry is the number of vertices with degree i. Try:`):
print(`DegSeq(GG(5),5);`):

elif nargs=1 and args[1]=Fc then
print(`Fc(c,x,y): sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2+c*x*y)`):

elif nargs=1 and args[1]=FEtoSeq then
print(`FEtoSeq(P,f,N): Given a polynomial P in f and a positive integer N, outputs the first N coefficients of`):
print(`the formal power seriers satisfying f(x)=x*P(f(x)). Try: `):
print(`FEtoSeq(1+f^2,f,10);`):

elif nargs=1 and args[1]=GG then
print(`GG(m,n,S): inputs positive integers m and n a subset S of [1,m]x[1,n] and outputs the graph in our data structure`):
print(`in the from [N,E,C] where N is the number of vertices (alias m*n-nops(S)) followed by the set of edges E ( a set of pairs in {1,...,N}, followed by the code C`):
print(`where C[n] is the actual grid point. Try:`):
print(`GG(6,6,{[6,1],[6,6]});`):

elif nargs=1 and args[1]=GGsp then
print(`GGsp(m,n,S,F): inputs positive integers m and n a subset S of [1,m]x[1,n], and a set of forbidden degrees F,and outputs `):
print(` a list of length 3 consisisting of`):
print(`(i) the set of vertices (as grid points [i,j]) using matrix notation [i,j] a subset of [1,m]x[1,n]`):
print(`(ii) the set of all edges`):
print(`(iii) the set of all spanning trees where the degrees are not allowed to be in F. Try:`):
print(`GGsp(5,5,{[1,4],[1,5],[5,5]},{2,4});`):



elif nargs=1 and args[1]=GP then
print(`GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)-3 in n,`):
print(`let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL.`):
print(`Try:`):
print(`GP([seq(i^3,i=1..10)],n);`):

elif nargs=1 and args[1]=InfoLT then
print(`InfoLT(S,K1,K2,n,N): inputs a list of positive integers S (that must contain 1),  positive integer K1, a large positive intger K2 and symbols n and N`):
print(`gives information about the sequence: Number of labeled trees where every vertex must have a degree that belongs to S`):
print(`The output is a list consisting of`):
print(`(i): The first K terms of the enumerating sequence`):
print(`(ii) the linear recurrence operator annihilating the sequence`):
print(`(iii) Just for one, the K2-th term`):
print(`(iv) the asympotics in the form [theta,mu] where the asymptotics is n!*mu^n*N^theta. Try:`):
print(`(v) the limiting distribution of the vertices according to their degrees (in S), followed by their respective variances.`):
print(`InfoLT([1,2,3],30,1000,n,N);`):


elif nargs=1 and args[1]=InfoLTeven then
print(`InfoLTeven(S,K1,C,K2,n,N): Like InfoLT (q.v.) but when the number of vertices has to be even to get a non-zero number. Try: `):
print(`InfoLTeven({1,3},30,10,1000,n,N);`):


elif nargs=1 and args[1]=InfoLTv then
print(`InfoLTv(S,K1,K2,n,N): A verbose form of InfoLT (q.v.). Try:`):
print(`InfoLTv([1,2,3],30,1000,n,N);`):

elif nargs=1 and args[1]=Kids then
print(`Kids(n,G,T,E,L): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves`):
print(`(it is implied by T and E, but it is convenient to have) outputs the set of all "children" such partial trees`):
print(`where one of the leaves gets children. Try:`):
print(`Kids(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4});`):


elif nargs=1 and args[1]=Kids1 then
print(`Kids1(n,G,T,E,L,a): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves`):
print(`(it is implied by T and E, but it is convenient to have) and a member a of L, outputs the set of all "children" such partial trees`):
print(`where a is no longer a leaf but has children. Try:`):
print(`Kids1(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4},2);`):


elif nargs=1 and args[1]=LimAveFor then
print(`LimAveFor(d): the formula, exp(-1)/(d-1)!, for the asympotic expercation of the r.v. "Number of verticecs in a random labeled tree on n vertices with degree d". Try:`):
print(`LimAveFor(d);`):

elif nargs=1 and args[1]=LimCorFor then
print(`LimCorFor(d1,d2): the formula for the limit of the correlation of the r.v.s "Number of verticecs in a random labeled tree on n vertices with degrees d1 and d2" divided by n. Try:`):
print(`LimCorFor(d1,d2);`):

elif nargs=1 and args[1]=LimCovFor then
print(`LimCovFor(d1,d2): the formula ,exp(-2)*(2*d1+2*d2-d1*d2-5)/((d1-1)!*(d2-1)!), for the limit of the covariance of the r.v.s "Number of verticecs in a random labeled tree on n vertices with degrees d1 and d2" divided by n. Try:`):
print(`LimCovFor(d1,d2);`):


elif nargs=1 and args[1]=LimMixedMoms then
print(`LimMixedMoms(d1,d2,L): The limit of the scaled momemts of the random variables "number of vertices with degree d1" and "number of vertices of degree d2" in labeled trees with n vertices, as n goes to infinity. Try:`):
print(`LimMixedMoms(1,2,4);`):


elif nargs=1 and args[1]=LimNor then
print(`LimNor(c,N): the list of lists of length N, let's call it L, such that its L[i][j] is the scaled momement of the joint normal distribution with`):
print(` pdf sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2-c*x*y). Try: `):
print(` LimNor(1/2,4); `):


elif nargs=1 and args[1]=LimVarFor then
print(`LimVarFor(d): the formula ,exp(-1)/(d-1)!- (d^2-4*d+5)/(d-1)!^2, for the limit of the variance of the r.v. "Number of verticecs in a random labeled tree on n vertices with degree d" divided by n. Try:`):
print(`LimVarFor(d);`):

elif nargs=1 and args[1]=LTFpaper then
print(`LTFpaper(M,k,K,n): inputs a positive integers M, and outputs the first K terms of the sequence enumerating the sequence`):
print(`number of labeled trees where no vertex has degree in the set S, for all subsets of {2,3,...,M} with at most k members.`):
print(`It also outputs the estimated asympotics of the form n^(n-2)*mu^n. It also outputs the list of the [S,mu] Try:`):
print(`LTFpaper(4,1,100,n);`):

elif nargs=1 and args[1]=LTmom then
print(`LTmom(P,n,a): inputs a list of positive integers P, a symbol n, and a list of positive integers a (of the same length of P) outputs the`):
print(`mixed-a-moments of the random variables (in the population of labelled trees with n vertices) "number of ocurrences of degree P[i] vertices.`):
print(`LTmom([1],n,[2]);`):

elif nargs=1 and args[1]=LTmom1 then
print(`LTmom1(d,n,a): The a-th moment (not about the mean) of the random variable "number of vertices of degree d" in the set of labeled trees on n vertices. Try:`):
print(`LTmom1(1,n,2);`):

elif nargs=1 and args[1]=LTmom1am then
print(`LTmom1am(d,n,a): The a-th moment (about the mean) of the random variable "number of vertices of degree d" in the set of labeled trees on n vertices. `):
print(`For example for the variance for the number of leaves, type:`):
print(`LTmom1am(1,n,2);`):


elif nargs=1 and args[1]=LTmom1amL then
print(`LTmom1amL(d,a): The limit of the a-th moment (about the mean) of the random variable "number of vertices of degree d"/n^(trunc(a/2)),  in the set of labeled trees on n vertices. `):
print(`For example for the limit of the variance for the number of leaves, type:`):
print(`LTmom1amL(1,2);`):


elif nargs=1 and args[1]=LTmom2 then
print(`LTmom2(d1,d2,n,a1,a2): The (a1,a2) mixed moments (not about the mean) of the random variables "number of vertices of degree d1"  abd "number of vertices with degree d2", in the set of labeled trees on n vertices. Try:`):
print(`LTmom2(1,2,n,2,2);`):


elif nargs=1 and args[1]=LTmom2am then
print(`LTmom2am(d1,d2,n,a1,a2): The (a1,a2) mixed moments, ABOUT  THE MEAN, of the random variables "number of vertices of degree d1"  abd "number of vertices with degree d2", in the set of labeled trees on n vertices. Try:`):
print(`LTmom2am(1,2,n,2,2);`):


elif nargs=1 and args[1]=LTmom2amL then
print(`LTmom2amL(d1,d2,a1,a2): The limit of LTmom2am(d1,d2,n,a1,a2)/n^((a1+a2)/2)) as n goes to infinity. Try:`):
print(`LTmom2amL(1,2,2,2);`):


elif nargs=1 and args[1]=LTpaper then
print(`LTpaper(M,K1,n): inputs a positive integer M, larger than 3, a positive integer K1, and a symbol n gives`):
print(`lots of information from InfoLT (q.v.) for all subsets of {1,...,M} with at least three elements and that contain 1. Try:`):
print(`LTpaper(3,30,n);`):


elif nargs=1 and args[1]=LTpaperPlain then
print(`LTpaperPlain(M,K1,n): inputs a positive integer M, larger than 3, a positive integer K1, and a symbol n gives`):
print(`lots of information from InfoLTplain (q.v.) for all subsets of {1,...,M} with at least three elements and that contain 1. Try:`):
print(`LTpaperPlain(3,30,n);`):

elif nargs=1 and args[1]=LTrec then
print(`LTrec(P,n,N): The linear recurrence operator annihilating the sequence enumerating labelled trees where all vertices have degrees that belong to P. Try:`):
print(`LTrec(P,n,N)`):

elif nargs=1 and args[1]=LTrecG then
print(`LTrecG(P,n,N,g): The linear recurrence operator annihilating the sequence weight-enumerating labelled trees where all vertices have degrees that belong to P, according to the weight Prod(g[degree(v)], over all vertices. Try:`):
print(`LTrecG(P,n,N,g);`):

elif nargs=1 and args[1]=LTseq then
print(`LTseq(P,N): The first N terms of the sequence enumerating labeled trees where all the vertex-degrees must belong to the set P. It uses the recurrence. Try:`):
print(`LTseq({1,2,3},10);`):


elif nargs=1 and args[1]=LTseqSlow then
print(`LTseqSlow(P,N): The first N terms of the sequence enumerating labeled trees where all the vertex-degrees must belong to the set P. It is much slower than LTseq(P,N), since it directly uses the functional equation.`):
print(`Nevertheless it is needed for the initial condtions of the recurrence. Try:`):
print(`LTseqSlow({1,2,3},10);`):


elif nargs=1 and args[1]=LTseqEven then
print(`LTseqEven(P,N): The first N terms of the sequence enumerating labeled trees, with 2*n vertices, in the case where the odd ones are all 0, where all the vertex-degrees must belong to the set P. Try:`):
print(`LTseqEven({1,3},10);`):

elif nargs=1 and args[1]=LTseqF then
print(`LTseqF(P,N): The first N terms of the sequence enumerating labeled trees where all the vertex-degrees must NOT belong to the set P`):
print(`Try:`):
print(`LTseqF({2},10);`):

elif nargs=1 and args[1]=LTseqFwe then
print(`LTseqFwe(P,N,g): The first N terms of the sequence weight-enumerating labeled trees according to the weight Prod(g[i]^(Number of vertices with degree i) in the list pf psoitive integers P. Try:`):
print(`LTseqFwe([1,2],10,g);`):

elif nargs=1 and args[1]=LTseqFwe1 then
print(`LTseqFwe1(P,n,g): The  weight-enumerator  labeled trees  with exactlty n vertices using the Lagrange inversion formula according to the weight Prod(g[i]^(Number of vertices with degree i) in the list pf psoitive integers P`):
print(`Try:`):
print(`LTseqFwe1([1,2],10,g);`):

elif nargs=1 and args[1]=LTseqWE then
print(`LTseqWE(P,N,g): The first N terms of the sequence weight-enumerating labeled trees where all the vertex-degrees must belong to the set P according to the weight : product of g[degree(v)] over all vertices v.`):
print(`Try:`):
print(`LTseqWE({1,2,3},30,g);`):



elif nargs=1 and args[1]=LTstat then
print(`LTstat(P,K1,r):Inputs a LIST,P,  of positive integers (including 1), of length at least 3,  and a  large integer K1, and a small integer r, outputs`):
print(`the ESTIMATED average/n, variance/n, and the scaled 3rd thoroug r-th moment of the random variables "number of vertices of degree P[i]"`):
print(`among labeled trees with n vertices, where all the vertices must be of degree in P`):
print(`Using the data for n=K1 and n=K1-1. Try:`):
print(`LTstat([1,2,3],1000,6);`):

elif nargs=1 and args[1]=LTstatJ1 then
print(`LTstatJ1(P,K1,r,deg1,deg2):Inputs a LIST of positive integers (including 1) and a  large integer K1, and a small integer r, `):
print(`and two members deg1,deg2, of P outputs scaled mixed moments for the random variables "number of vertices with degree deg1`):
print(`and "number of vertices of degree deg2" up to the r-th mixed moment. Try:`):
print(`LTstatJ1([1,2,3],1000,2,1,2);`):

elif nargs=1 and args[1]=MixMomNorF then
print(`MixMomNorF(c,m,n): The mixed-moment of the bivairate normal distribution with correlation c. Done the fast way, using recurrences obtained via the multivariable Almkvist-Zeilberger. Try:`):
print(`MixNorF(c,4,6);`)


elif nargs=1 and args[1]=ProveAsyNor then
print(`ProveAsyNor(d1,r): Rigorously proves the asymptotic normality of the random varialbe "number of vertices with degree d" using the method of moments by checking that the scaled ven moments are the same as those of the normal distribution`):
print(`i.e. 1,3,15,105, ... to the 2r-th moment. Try:`):
print(`ProveAsyNor(1,5);`):

elif nargs=1 and args[1]=ProveJointAsyNor then
print(`ProveJointAsyNor(d1d2,L): Rigorously proves the joint asymptotic (but, not independent) normality of the random varialbes `):
print(`"number of vertices with degree d1"  and "number of vertices with degree d2" in labeled trees`):
print(`using the method of moments by checking that the scaled mixed moments are the same as those of the scaled mixed `):
print(`moments with the same correlation. Try`):
print(`ProveJointAsyNor(1,2,4);`):

elif nargs=1 and args[1]=RG then
print(`RG(n,K): A random graph with n vertices and K edges. Try:`):
print(`RG(5,12);`):


elif nargs=1 and args[1]=ScaledMixedMoms then
print(`ScaledMixedMoms(f,x,y,N): Given a probability generating function`):
print(`f(x,y) (a polynomial, rational function and possibly beyond)`):
print(`returns a double L such that`):
print(`L[1][1] is the correlation  (in decimals) and L[i][j] is the scaled (i,j) moment`):
print(`For example, try:`):
print(`ScaledMixedMoments(((x+y+1)^100,x,y,4);`):


elif nargs=1 and args[1]=StatPaper then
print(`StatPaper(MaxDeg,MaxMom1, MaxMom2): Inputs positive integers MaxDeg, MaxMom1, and MaxMom2 and outputs a paper about the r.v. "number of vertices of degree d"`):
print(`with explicit expressions for the moments for d between 1 amd MaxMom1 and for the mixed moments up to  MaxMom2. Try:`):
print(`StatPaper(4,4,4);`):

elif nargs=1 and args[1]=SubsFPS then
print(`SubsFPS(L,P,f): Given a sequence of numbers L, and a polynomial P in f, let f(x)=Sum(L[i]*x^i,i=0..N). Outputs the list consisting of the first N coefficients`):
print(`of x*P(f(x)). Try:`):
print(`SubsFPS([seq(binomial(2*i,i)/(i+1),i=1..10)],1+f^2,f);`):

elif nargs=1 and args[1]=STmtt then
print(`STmtt(n,G): Given a simple graph G on n vertices outputs the set of spanning trees. It uses the Matrix Tree Theorem and computer algebra.`):
print(`Try: `):
print(`STmtt(5,RG(5,10));`):


elif nargs=1 and args[1]=ST then
print(`ST(n,G): Given a simple graph G on n vertices outputs the set of spanning trees. It uses dynamical programming.`):
print(`Try: `):
print(`ST(3,RG(3,3));`):


elif nargs=1 and args[1]=STf then
print(`STf(n,G,F): Given a simple graph G on n vertices and a set of forbidden degrees, F, outputs the set of spanning trees where no vertex can have a degree in F. It uses dynamical programming.`):
print(`Try: `):
print(`STf(3,RG(3,3),{2});`):


elif nargs=1 and args[1]=STnu then
print(`STnu(n,G): Given a simple graph G on n vertices outputs the NUMBER of spanning trees. It uses the Matrix Tree Theorem`):
print(`Try: `):
print(`STnu(5,RG(5,10));`):

elif nargs=1 and args[1]=STd then
print(`STd(n,G,d): Given a simple graph G on n vertices and a degree sequence of length n, outputs the set of spanning trees with the degree sequence. Try:`):
print(`STd(3,{{1,2},{1,3},{2,3}},[1,1,2]);`):

elif nargs=1 and args[1]=STmttF then
print(`STmttF(n,G,F): Same as STf(n,G,F) but using the Matrix Tree Theorem and symbolic computations.`):
print(` Try`):
print(`STmttF(5,RG(5,10),{2});`):


elif nargs=1 and args[1]=STmttP then
print(`STmttP(n,G,P): Same as STp(n,G,P) but using the Matrix Tree Theorem and symbolic computations.`):
print(` Try`):
print(`STmttP(5,RG(5,10),{1,3});`):


elif nargs=1 and args[1]=STfNu then
print(`STfNu(n,G,F): Given a simple graph G on n vertices and a set of positive integers F, outputs  the NUMBER of spanning trees of G`):
print(`where none of the degrees of the vertices belongs to F. Try`):
print(`STfNu(5,RG(5,10),{2,4});`):


elif nargs=1 and args[1]=STp then
print(`STp(n,G,P): Given a simple graph G on n vertices and a set of allowed degrees, P, outputs the set of spanning trees where every vertex has a degree in P. It uses dynamical programming.`):
print(`Try: `):
print(`STp(5,CG(5),{2});`):

elif nargs=1 and args[1]=T1rn then
print(`T1rn(r,n): the number of labeled trees on n vertices where none of the vertices as degree r. Try:`):
print(`T1rn(2,20);`):


elif nargs=1 and args[1]=T2rn then
print(`T2rn(r1,r2,n): the number of labeled trees on n vertices where none of the vertices has a degree in the set {r1,r2}. Try:`):
print(`T2rn(2,3,20);`):


elif nargs=1 and args[1]=T3rn then
print(`T3rn(r1,r2,r3,n): the number of labeled trees on n vertices where none of the vertices has a degree in the set {r1,r2,r3}. Try:`):
print(`T3rn(2,3,4,20);`):

elif nargs=1 and args[1]=Wt then
print(`Wt(G,n,g): The weight of the graph G according to product(vertices) g[i]^degree[i]. Try`):
print(`Wt({{1,2},{1,3}},3,g);`):

elif nargs=1 and args[1]=Zinn then
print(` Zinn(L): inputs a sequence L,  uses Zinn-Justin's method to approximate mu and theta such that L[n] is asymptotically C*n^theta*mu^n for some constant C. Try`):
print(`Zinn([seq(2^i*i,i=1..100)]);`):

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

fi:


end:




with(combinat):
with(linalg):


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:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
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:


##Start from GuessPol

#GP1(L,n,d): inputs a list L , a variable name n, and a positive integer d, guesses a polynomial of degree <=d in n,
#let's call it P(n) such that P(i)=L[i]. The length of L must be at least d+3. If it fails it returns FAIL.
#Try:
#GP1([seq(i^3,i=1..10)],n,3);
GP1:=proc(L,n,d) local a,i,var,eq,P:

#Any d+1 data points can fit a polynomial of degree d. We want at least 2 "confirmation" points
if nops(L)<=d+2 then
 print(`The list must be of length at least`, d+3):
 RETURN(FAIL):
fi:

#We define a generic polynomial P of n, of degree d, with UNDETERMINED coefficients. Our goal is to determine them

P:=add(a[i]*n^i,i=0..d):

#This is the set of unknowns with d+1 unknowns
var:={seq(a[i],i=0..d)}:

#We set-up a system with d+3 equations fitting the polynomial to the data
eq:={seq(subs(n=i,P)-L[i],i=1..d+3)}:

#Maple tries to solve the system, we rename var the solution
var:=solve(eq,var):

#If var=NULL then there is no solution, in other words, no polynomial of degree d fits the data

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

#Now we plug the solution var into the generic polynomial P, making it an actual polynomial (that fits the data)
P:=subs(var,P):

#If nops(L) is larger than d+3, we check that the polynomial fits the remaining data
for i from d+4 to nops(L) do
 if expand(subs(n=i,P))-L[i]<>0 then
  print(P, `does not agree with n=`, i):
  RETURN(FAIL):
 fi:
od:

#If it agrees we return the polynomial
P:

end:


  

#GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)-3 in n,
#let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL.
#Try:
#GP([seq(i^3,i=1..10)],n);

GP:=proc(L,n) local d,P:

for d from 0 to nops(L)-3 do
#We start with d=0 and keep trying with higher and higher degree
 P:=GP1(L,n,d):
 if P<>FAIL then
  RETURN(P):
 fi:
od:
#If there is no such polynomial of degree <=nopd(L)-3 it returns FAIL
FAIL:
end:

##end from GuessPol
#FROM Findrec

ezraFindrec:=proc()
if args=NULL then

print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

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

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

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


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


fi:
 
end:

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

ope:=subs(var1,ope):

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

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

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

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


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

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

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


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

ope:=subs(var1,ope):

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


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

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

end:
 


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


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L,n1,i:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
if {seq(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..degree(ope,N)),n1=1..nops(f)-degree(ope,N))}<>{0} then
 RETURN(FAIL):
fi:
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





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

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

gu:=Ini:

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

gu:

end:
 
#end Findrec


with(linalg):

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


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


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

A1:=matrix(D1,D1):

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

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

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

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

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

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

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

true:

end:


######end from Findrec

#CG(n): A random graph with n vertices and K edges
CG:=proc(n) local i,j:{seq(seq({i,j},j=i+1..n),i=1..n)}:end:

#RG(n,K): A random graph with n vertices and K edges
RG:=proc(n,K) local i,j:

if not (type(n,posint) and type(K,posint) and K<=binomial(n,2)) then
 print(`bad input`):
 RETURN(FAIL):
fi:

randcomb({seq(seq({i,j},j=i+1..n),i=1..n)},K):


end:


#STmtt(n,G): Given a simple graph G on n vertices outputs the set of spanning trees. It uses the Matrix Tree theorem

STmtt:=proc(n,G) local i,j,M,e,lu,a,i1,j1:

if not (type(n,integer) and n>0 and type(G,set)) then
 print(`bad input`):
 RETURN(FAIL):
fi:

for i from 1 to n-1 do
for j from 1 to n-1 do
M[i,j]:=0:
od:
od:

for e in G do
 M[e[1],e[2]]:=M[e[1],e[2]]-a[e]:
 M[e[2],e[1]]:=M[e[2],e[1]]-a[e]:
 M[e[1],e[1]]:= M[e[1],e[1]]+a[e]:
 M[e[2],e[2]]:= M[e[2],e[2]]+a[e]:
od:

lu:=expand(det([seq([seq(M[i,j],j=1..n-1)],i=1..n-1)])):

if lu=0 then
 RETURN({}):
fi:

if type(lu,`+`) then
{seq({seq(op(1,op(i1,op(j1,lu))),i1=1..nops(op(1,lu)))},j1=1..nops(lu))}:

elif type(lu,`*`) then

{{seq(op(1,op(i1,lu)),i1=1..nops(lu))}}:

else
FAIL:
fi:


end:


#STd(n,G,d): Given a simple graph G on n vertices and a degree sequence of length n, outputs the set of spanning trees with the degree sequence. Try:
#STd(3,{{1,2},{1,3},{2,3}},[1,1,2]);
STd:=proc(n,G,d) local i,j,M,e,lu,a,x,i1,j1:

if not (type(n,integer) and n>0 and type(G,set) and nops(d)=n) and convert(d,`+`)=2*(n-1) then
 print(`bad input`):
 RETURN(FAIL):
fi:


for i from 1 to n-1 do
for j from 1 to n-1 do
M[i,j]:=0:
od:
od:

for e in G do
 M[e[1],e[2]]:=M[e[1],e[2]]-a[e]*x[e[1]]*x[e[2]]:
 M[e[2],e[1]]:=M[e[2],e[1]]-a[e]*x[e[1]]*x[e[2]]:
 M[e[1],e[1]]:= M[e[1],e[1]]+a[e]*x[e[1]]*x[e[2]]:
 M[e[2],e[2]]:= M[e[2],e[2]]+a[e]*x[e[1]]*x[e[2]]:
od:

lu:=expand(det([seq([seq(M[i,j],j=1..n-1)],i=1..n-1)])):




for i from 1 to n do
 lu:=coeff(lu,x[i],d[i]):
od:


if lu=0 then
 RETURN({}):
fi:

if type(lu,`+`) then
{seq({seq(op(1,op(i1,op(j1,lu))),i1=1..nops(op(1,lu)))},j1=1..nops(lu))}:

elif type(lu,`*`) then

{{seq(op(1,op(i1,lu)),i1=1..nops(lu))}}:

else
FAIL:
fi:


end:


#STmttF(n,G,F): Given a simple graph G on n vertices and a set of positive integers F, outputs all the spanning trees of G
#where none of the degrees of the vertices belongs to F. Try
#SFmttF(5,RG(5,10),{2,4});

STmttF:=proc(n,G,F) local i,j,M,e,lu,a,x,f,i1,j1:

if not (type(n,integer) and n>0 and type(G,set) and type(F,set) and (F={} or {seq(type(f,posint),f in F)}={true})) then
 print(`bad input`):
 RETURN(FAIL):
fi:


for i from 1 to n-1 do
for j from 1 to n-1 do
M[i,j]:=0:
od:
od:

for e in G do
 M[e[1],e[2]]:=M[e[1],e[2]]-a[e]*x[e[1]]*x[e[2]]:
 M[e[2],e[1]]:=M[e[2],e[1]]-a[e]*x[e[1]]*x[e[2]]:
 M[e[1],e[1]]:= M[e[1],e[1]]+a[e]*x[e[1]]*x[e[2]]:
 M[e[2],e[2]]:= M[e[2],e[2]]+a[e]*x[e[1]]*x[e[2]]:
od:

lu:=expand(det([seq([seq(M[i,j],j=1..n-1)],i=1..n-1)])):




for i from 1 to n do
 for f in F do
  lu:=expand(lu-coeff(lu,x[i],f)*x[i]^f):
 od:
od:

lu:=subs({seq(x[i]=1,i=1..n)},lu):


if lu=0 then
 RETURN({}):
fi:

if type(lu,`+`) then
{seq({seq(op(1,op(i1,op(j1,lu))),i1=1..nops(op(1,lu)))},j1=1..nops(lu))}:

elif type(lu,`*`) then

{{seq(op(1,op(i1,lu)),i1=1..nops(lu))}}:

else
FAIL:
fi:


end:




#STmttP(n,G,P): Given a simple graph G on n vertices and a set of positive integers P, outputs all the spanning trees of G
#where every degree  belongs to P. Try
#SFmttP(5,RG(5,10),{1,3});

STmttP:=proc(n,G,P) local i,j,M,e,lu,a,x,p,i1,j1,mu,lu1,f:

if not (type(n,integer) and n>0 and type(G,set) and type(P,set) and (member(1,P) and  {seq(type(f,posint),f in P)}={true})) then
 print(`bad input`):
 RETURN(FAIL):
fi:


for i from 1 to n-1 do
for j from 1 to n-1 do
M[i,j]:=0:
od:
od:

for e in G do
 M[e[1],e[2]]:=M[e[1],e[2]]-a[e]*x[e[1]]*x[e[2]]:
 M[e[2],e[1]]:=M[e[2],e[1]]-a[e]*x[e[1]]*x[e[2]]:
 M[e[1],e[1]]:= M[e[1],e[1]]+a[e]*x[e[1]]*x[e[2]]:
 M[e[2],e[2]]:= M[e[2],e[2]]+a[e]*x[e[1]]*x[e[2]]:
od:

lu:=expand(det([seq([seq(M[i,j],j=1..n-1)],i=1..n-1)])):


mu:=0:

for lu1 in lu do

  if {seq(degree(lu1,x[i1]),i1=1..n)} minus P={} then
   mu:=mu+lu1:
  fi:
od:



mu:=subs({seq(x[i]=1,i=1..n)},mu):


if mu=0 then
 RETURN({}):
fi:

if type(mu,`+`) then
{seq({seq(op(1,op(i1,op(j1,mu))),i1=1..nops(op(1,mu)))},j1=1..nops(mu))}:

elif type(lu,`*`) then

{{seq(op(1,op(i1,mu)),i1=1..nops(mu))}}:

else
FAIL:
fi:


end:



#GG(m,n,S): inputs positive integers m and n a subset S of [1,m]x[1,n] and outputs the graph in our data structure
#in the from [N,E,C] where N is the number of vertices (alias m*n-nops(F)) followed by the set of edges E ( a set of pairs in {1,...,N}, followed by the code C
#where C[n] is the actual grid point. Try:
#GG(6,6,{[6,1],[6,6]});
GG:=proc(m,n,S) local V,E,i,j,v,Nei,nei:
V:={seq(seq([i,j],j=1..n),i=1..m)} minus S:
E:={}:

for v in V do

 Nei:={v+[1,0],v-[1,0],v+[0,1],v-[0,1]} intersect V:
 E:=E union {seq({v,nei},nei in Nei)}:
od:

V,E:

V:=convert(V,list):
E:=subs({seq(V[i]=i,i=1..nops(V))},E):
[nops(V),E,V]:

end:




#STfNu(n,G,F): Given a simple graph G on n vertices and a set of positive integers F, outputs  the NUMBER of spanning trees of G
#where none of the degrees of the vertices belongs to F. Try
#SFfNu(5,RG(5,10),{2,4});

STfNu:=proc(n,G,F) local i,j,M,e,lu,x,f:

if not (type(n,integer) and n>0 and type(G,set) and type(F,set) and (F={} or {seq(type(f,posint),f in F)}={true})) then
 print(`bad input`):
 RETURN(FAIL):
fi:


for i from 1 to n-1 do
for j from 1 to n-1 do
M[i,j]:=0:
od:
od:

for e in G do
 M[e[1],e[2]]:=M[e[1],e[2]]-x[e[1]]*x[e[2]]:
 M[e[2],e[1]]:=M[e[2],e[1]]-x[e[1]]*x[e[2]]:
 M[e[1],e[1]]:= M[e[1],e[1]]+x[e[1]]*x[e[2]]:
 M[e[2],e[2]]:= M[e[2],e[2]]+x[e[1]]*x[e[2]]:
od:

lu:=expand(det([seq([seq(M[i,j],j=1..n-1)],i=1..n-1)])):




for i from 1 to n do
 for f in F do
  lu:=expand(lu-coeff(lu,x[i],f)*x[i]^f):
 od:
od:

subs({seq(x[i]=1,i=1..n)},lu):



end:


#STnu(n,G): Given a simple graph G on n vertices. Finds the number of spanning trees
#STnu(5,RG(5,10),{2,4});
#STnu(34,GG(6,6,{[6,1],[6,6]})[2]):
STnu:=proc(n,G) local i,j,M,e:

if not (type(n,integer) and n>0 and type(G,set)) then
 print(`bad input`):
 RETURN(FAIL):
fi:


for i from 1 to n-1 do
for j from 1 to n-1 do
M[i,j]:=0:
od:
od:

for e in G do
 M[e[1],e[2]]:=M[e[1],e[2]]-1:
 M[e[2],e[1]]:=M[e[2],e[1]]-1:
 M[e[1],e[1]]:= M[e[1],e[1]]+1:
 M[e[2],e[2]]:= M[e[2],e[2]]+1:
od:

det([seq([seq(M[i,j],j=1..n-1)],i=1..n-1)]);


end:


#Kids1(n,G,T,E,L,a): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves
#(it is implied by T and E, but it is convenient to have) and a member a of L, outputs the set of all "children" such partial trees
#where a is no longer a leaf but has children. Try:
#Kids1(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4},2);
Kids1:=proc(n,G,T,E,L,a) local AvN,v,s,gu,T1,E1,L1,i:

if not member(a,L) then
 RETURN(FAIL):
fi:

AvN:={}:

for v in {seq(i,i=1..n)} minus T do
 if member({a,v}, G minus E) then
  AvN:=AvN union {v}:
 fi:
od:

if AvN={} then
 RETURN({}):
fi:

AvN:=powerset(AvN) minus {{}}:

gu:={}:

for s in AvN do
 T1:=T union s:
 E1:=E union {seq({a,v}, v in s)}:
 L1:=(L minus {a}) union s:
 gu:=gu union {[n,G,T1,E1,L1]}:
od:

gu:

end:




 


#Kids(n,G,T,E,L): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves
#(it is implied by T and E, but it is convenient to have) outputs the set of all "children" such partial trees
#where one of the leaves gets children. Try:
#Kids(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4});
Kids:=proc(n,G,T,E,L) local a:
{seq(op(Kids1(n,G,T,E,L,a)), a in L)}:
end:




#ST(n,G): Given a simple graph G on n vertices outputs the set of spanning trees, the fast way

ST:=proc(n,G) local gu,cu,cuNew,cu1:

cu:={[n,G,{1},{},{1}]}:

gu:={}:

cu:=Kids(op(cu[1])):


 for cu1 in cu do
  if nops(cu1[4])=n-1 then
   gu:=gu union {cu1[4]}:
  fi:
od:

while cu<>{} do
cuNew:={}:

for cu1 in cu do
  cuNew:=cuNew union Kids(op(cu1)):
od:

cu:=cuNew:

 for cu1 in cu do
  if nops(cu1[4])=n-1 then
   gu:=gu union {cu1[4]}:
  fi:
od:

od:
gu:

end:





###start with forbidden degrees



#Kids1F(n,G,T,E,L,a,F): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves
#(it is implied by T and E, but it is convenient to have) and a member a of L, and a set of forbidden number of children, outputs the set of all "children" such partial trees
#where a is no longer a leaf but has children. Try:
#Kids1F(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4},2,{1,3});
Kids1F:=proc(n,G,T,E,L,a,F) local AvN,v,s,gu,T1,E1,L1,i:

if not member(a,L) then
 RETURN(FAIL):
fi:

AvN:={}:

for v in {seq(i,i=1..n)} minus T do
 if member({a,v}, G minus E) then
  AvN:=AvN union {v}:
 fi:
od:

if AvN={} then
 RETURN({}):
fi:

AvN:=powerset(AvN) minus {{}}:

gu:={}:

for s in AvN do

 if not member(nops(s),F) then

 T1:=T union s:
 E1:=E union {seq({a,v}, v in s)}:
 L1:=(L minus {a}) union s:
 gu:=gu union {[n,G,T1,E1,L1]}:
 fi:
od:

gu:

end:


#KidsF(n,G,T,E,L,F): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves
#(it is implied by T and E, but it is convenient to have) and a forbidden number of children, outputs the set of all "children" such partial trees
#where one of the leaves gets children. Try:
#KidsF(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4},{2,4});
KidsF:=proc(n,G,T,E,L,F) local a:
{seq(op(Kids1F(n,G,T,E,L,a,F)), a in L)}:
end:


#STf(n,G,F): Given a simple graph G on n vertices, and a set of forbidden degrees, outputs the set of spanning trees, the fast way

STf:=proc(n,G,F) local gu,cu,cuNew,cu1,F1,i1:

if member(1,F) then
 RETURN({}):
fi:

cu:={[n,G,{1},{},{1}]}:

gu:={}:

cu:=KidsF(op(cu[1]),F):


F1:={seq(i1-1,i1 in F)}:

 for cu1 in cu do
  if nops(cu1[4])=n-1 then
   gu:=gu union {cu1[4]}:
  fi:
od:

while cu<>{} do
cuNew:={}:

for cu1 in cu do
  cuNew:=cuNew union KidsF(op(cu1),F1):
od:

cu:=cuNew:

 for cu1 in cu do
  if nops(cu1[4])=n-1 then
   gu:=gu union {cu1[4]}:
  fi:
od:

od:
gu:

end:



#GGsp(m,n,S,F): inputs positive integers m and n a subset S of [1,m]x[1,n], and a set of forbidden degrees F,and outputs 
# a list of length 3 consisisting of
#(i) the set of vertices (as grid points [i,j]) using matrix notation [i,j] a subset of [1,m]x[1,n]
#(ii) the set of all edges
#(iii) the set of all spanning trees where the degrees are not allowed to be in F. Try:
#GGsp(5,5,{[1,4],[1,5],[5,5]},{2,4});
GGsp:=proc(m,n,S,F)  local gu,V,G,i,lu:

gu:=GG(m,n,S):

V:=gu[3]:

G:=gu[2]:

lu:=STf(gu[1],G,F):

G:=subs({seq(i=V[i],i=1..nops(V))},G):

lu:=subs({seq(i=V[i],i=1..nops(V))},lu):


[{op(V)},G,lu]:

end:


####start spanning tree where only a finite set is permitted




#Kids1P(n,G,T,E,L,a,P): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves
#(it is implied by T and E, but it is convenient to have) and a member a of L, and a set of allowed number of children, outputs the set of all "children" such partial trees
#where a is no longer a leaf but has children. Try:
#Kids1P(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4},2,{2});
Kids1P:=proc(n,G,T,E,L,a,P) local AvN,v,s,gu,T1,E1,L1,i:

if not member(a,L) then
 RETURN(FAIL):
fi:

AvN:={}:

for v in {seq(i,i=1..n)} minus T do
 if member({a,v}, G minus E) then
  AvN:=AvN union {v}:
 fi:
od:

if AvN={} then
 RETURN({}):
fi:

AvN:=powerset(AvN) minus {{}}:

gu:={}:

for s in AvN do

 if  member(nops(s),P) then

 T1:=T union s:
 E1:=E union {seq({a,v}, v in s)}:
 L1:=(L minus {a}) union s:
 gu:=gu union {[n,G,T1,E1,L1]}:
 fi:
od:

gu:

end:


#KidsP(n,G,T,E,L,P): given a pos. integer n, a graph G, a partial spanning tree with a set of vertices T and a set of edges E, and a subset of T, the current leaves
#(it is implied by T and E, but it is convenient to have) and a permitted number of children, outputs the set of all "children" such partial trees
#where one of the leaves gets children. Try:
#KidsP(4,{{1,2},{1,3},{1,4},{2,3},{2,4},{3,4}},{1,2,4},{{1,2},{1,4}},{2,4},{3});
KidsP:=proc(n,G,T,E,L,P) local a:
{seq(op(Kids1P(n,G,T,E,L,a,P)), a in L)}:
end:

#STp(n,G,P): Given a pos. integer n and a graph G on {1,...,n} and a set of positive integers P larger than 1 outputs the set of
#spanning trees where every vertex is either a leaf or has a degree that belongs to P. Try:
#STp(5,CG(5),{2});
STp:=proc(n,G,P) local gu,cu,cuNew,cu1,P1,i1:

if not member(1,P) then
 RETURN(FAIL):
fi:


cu:={[n,G,{1},{},{1}]}:

gu:={}:

cu:=KidsP(op(cu[1]),P):


P1:={seq(i1-1,i1 in P)}:

 for cu1 in cu do
  if nops(cu1[4])=n-1 then
   gu:=gu union {cu1[4]}:
  fi:
od:

while cu<>{} do
cuNew:={}:

for cu1 in cu do
  cuNew:=cuNew union KidsP(op(cu1),P1):
od:

cu:=cuNew:

 for cu1 in cu do
  if nops(cu1[4])=n-1 then
   gu:=gu union {cu1[4]}:
  fi:
od:

od:
gu:

end:


#FEtoSeq(P,f,N): Given a polynomial P in f and a positive integer N, outputs the first N coefficients of
#the formal power seriers satisfying f(x)=x*P(f(x)). Try: 
#FEtoSeq(1+f^2,f,10);
FEtoSeq:=proc(P,f,N) local x,gu1,gu2,gu3,i:
gu1:=x:
gu2:=subs(f=gu1,x*P):
gu2:=taylor(gu2,x=0,N+1):
gu2:=expand(add(coeff(gu2,x,i)*x^i,i=1..N)):
while gu1<>gu2 do
gu3:=subs(f=gu2,x*P):
gu3:=taylor(gu3,x=0,N+1):
gu3:=expand(add(coeff(gu3,x,i)*x^i,i=1..N)):
gu1:=gu2:
gu2:=gu3:
od:
expand([seq(coeff(gu2,x,i),i=1..N)]):
end:


#SubsFPS(L,P,f): Given a sequence of numbers L, and a polynomial P in f, let f(x)=Sum(L[i]*x^i,i=0..N). Outputs the list consisting of the first N coefficients
#of x*P(f(x)). Try:
#SubsFPS([seq(binomial(2*i,i)/(i+1),i=1..10)],1+f^2,f);

SubsFPS:=proc(L,P,f)  local gu,x,i:
gu:=add(L[i]*x^i,i=1..nops(L)):
gu:=taylor(subs(f=gu,x*P),x=0,nops(L)+3):
expand([seq(coeff(gu,x,i),i=1..nops(L))]):
end:





####EVEN parts
#LTseqEven(P,N): The first N terms of the sequence enumerating labeled trees with 2n vertices, where all the vertex-degrees must belong to the set P
LTseqEven:=proc(P,N) local gu,i:

gu:=LTseq(P,2*N):
if gu=FAIL then
 RETURN(FAIL):
fi:

if {seq(gu[2*i-1],i=1..N)}<>{0} then
 print(`Not  only even `):
 RETURN(FAIL):
fi:
[seq(gu[2*i],i=1..N)]:
end:


#InfoLTeven(S,K1,K2,C,n,N): inputs a set of positive integers S (that must contain 1), a positive integers K1, and K2 (K2 much larger than K1),and symbols n and N
#gives information about the sequence: Number of labeled trees with 2*n vertices where every vertex must have a degree that belongs to S
#The output is a list consisting of
#(i): The first K1 terms of the enumerating sequence
#(ii) the linear recurrence operator annihilating the sequence if its complexity is <=C
#(iii) Just for fun, the K2-th term
#InfoLTeven({1,2,3},30,10,1000,2,10,n,N);
InfoLTeven:=proc(S,K1,C,K2,n,N) local gu,ope,gadol,K,gu1,i,ope1,lu,j:

gu1:=LTseqEven(S,K1):


K:=max(seq((2+i)*(1+C-i)+4,i=0..C))+10:

gu:=LTseqEven(S,K):




ope:=Findrec(gu,n,N,C):

if ope=FAIL then
 RETURN([gu1]):
fi:




gadol:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],K2)[K2]:

ope1:=add(coeff(ope,N,i)*mul(n+j,j=1..i)*N^i,i=0..degree(ope,N)):


lu:=As(ope1,n,N):

if lu=FAIL then
 [gu1,ope,gadol]:
else
 [gu1,ope,gadol,lu,evalf(lu)]:
fi:

end:

####End EVEN parts

#LTseqSlow(P,N): The first N terms of the sequence enumerating labeled trees where all the vertex-degrees must belong to the set P
LTseqSlow:=proc(P,N) local i,L,f:

if not (type(N,posint) and type(P,set) and member(1,P) and {seq(type(i,integer), i in P)}={true}) then
print(`bad input`):
 RETURN(FAIL):
fi:

L:=FEtoSeq(add(f^(i-1)/(i-1)!, i in P),f,N):
L:=SubsFPS(L,add(f^i/i!,i in P),f):
[seq(L[i]*(i-1)!,i=1..N)]:


end:




#LTseq(P,N): Same as LTseq(P,N) (q.v.) but must faster, using the recurrence. Try
#LTseq({1,2,3},1000)[1000];
LTseq:=proc(P,N) local i,L,f,ope,INI,n,Sn:

if not (type(P,set) and type(N,posint)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

ope:=LTrec(P,n,Sn):

INI:=LTseqSlow(P,degree(ope,Sn)):

if N<nops(INI) then
[op(1..N,INI)]:
else
SeqFromRec(ope,n,Sn,INI,N):
fi:
end:


#LTseqF(P,N): The first N terms of the sequence enumerating labeled trees where all the vertex-degrees must belong to the set P
#Try:
#LTseqF({2},10);
LTseqF:=proc(P,N) local i,L,f,P1,n1:


if not (type(N,posint) and type(P,set) and {seq(type(i,integer), i in P)}={true}) then
print(`bad input`):
 RETURN(FAIL):
fi:

if member(1,P) then
  RETURN(FAIL):
fi:

P1:=sort([op(P)]):
if nops(P1)=1 then
 RETURN([seq(T1rn(P1[1],n1),n1=1..N)]):
elif nops(P1)=2 then
 RETURN([seq(T2rn(P1[1],P1[2],n1),n1=1..N)]):
elif nops(P1)=3 then
 RETURN([seq(T3rn(P1[1],P1[2],P1[3],n1),n1=1..N)]):
else
L:=FEtoSeq(exp(f)- add(f^(i-1)/(i-1)!, i in P),f,N):
L:=SubsFPS(L,exp(f)-add(f^i/i!,i in P),f):
RETURN([0,seq(L[i]*(i-1)!,i=2..N)]):
fi:

end:





#As(ope,n,N): Outputs the connective constant and criticacl exponent
As:=proc(ope,n,N) local alpha,ope1,t,i,lu,a,f,k1,g,aluf,si,i1:

alpha:=[solve(lcoeff(numer(ope),n),N)]:

if alpha=[] then
 RETURN(FAIL):
fi:

aluf:=alpha[1]:
si:=abs(evalf(aluf)):

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

alpha:=aluf:


if abs(coeff(evalf(alpha),I,1))>10^(-7) or evalf(alpha)<0 then
 RETURN(FAIL):
fi:

ope1:=subs(N=N*alpha,ope):

lu:=numer(add(subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^t,i=0..degree(ope1,N))):

lu:=taylor(lu,x=0,31);

lu:=add(simplify(expand(coeff(lu,x,i)))*x^i,i=0..30):

for i1 from 1 to 29 while coeff(lu,x,i1)=0 do od:

if i1=30 then
 RETURN([alpha]):
else

t:=solve(coeff(lu,x,i1),t);

RETURN([alpha,t]):
fi:

end:



 
pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:
 
goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:
 
duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1:
KAMA:=10:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:


AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=20:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:
 




#LTrec(P,n,N): The linear recurrence operator annihilating the sequence enumerating labelled trees where all vertices have degrees that belong to P. Try:
#LTrec(P,n,N)
LTrec:=proc(P,n,N) local z,G,i,ope,c:
option remember:
G:=add(z^(i-1)/(i-1)!, i in P):
ope:=AZd(G^n/z^(n-1)*(n-2)!,z,n,N)[1]:
c:=lcoeff(ope,N):
add(factor(normal(coeff(ope,N,i)/c))*N^i,i=0..degree(ope,N)):
end:

#InfoLT(S,K1,K2,n,N): inputs a list of positive integers, S (that must contain 1), a positive integer K, and symbols n and N
#gives information about the sequence: Number of labeled trees where every vertex must have a degree that belongs to S
#The output is a list consisting of
#(i): The first K1 terms of the enumerating sequence
#(ii) the linear recurrence operator annihilating the sequence
#(iii) Just for fun, the K2-th term
#(iv): The  limiting distribution of the number of vertices with the degrees in S, followed by the respective variance (divided by n)
#InfoLT([1,2,3],30,10,1000,2,10,n,N);
InfoLT:=proc(S,K1,K2,n,N) local gu,ope,gadol,K,gu1,i,ope1,lu,j,ku,i1:



ope:=LTrec(convert(S,set),n,N):

gu1:=LTseq(convert(S,set),max(K1,degree(ope,N))):

gadol:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu1)],K2)[K2]:

ope1:=add(coeff(ope,N,i)*mul(n+j,j=1..i)*N^i,i=0..degree(ope,N)):


ku:=LTstat(S,300,2):

lu:=As(ope1,n,N):

if lu=FAIL and ku=FAIL then
 RETURN([gu1,ope,gadol,FAIL,FAIL,FAIL,FAIL ]):
elif lu<>FAIL and ku=FAIL then
RETURN([gu1,ope,gadol,lu,evalf(lu), FAIL,FAIL ]):
elif lu=FAIL and ku<>FAIL then
RETURN([gu1,ope,gadol,FAIL,FAIL, [seq(ku[i1][1],i1=1..nops(ku))] ,[seq(ku[i1][2],i1=1..nops(ku))] ]):
else 
RETURN([gu1,ope,gadol,lu,evalf(lu), [seq(ku[i1][1],i1=1..nops(ku))] ,[seq(ku[i1][2],i1=1..nops(ku))] ]):
fi:



end:


#InfoLTplain(S,K1,K2,n,N): inputs a list of positive integers, S (that must contain 1), a positive integer K, and symbols n and N
#gives information about the sequence: Number of labeled trees where every vertex must have a degree that belongs to S
#The output is a list consisting of
#(i): The first K1 terms of the enumerating sequence
#(ii) the linear recurrence operator annihilating the sequence
#(iii) Just for fun, the K2-th term
#InfoLTplain([1,2,3],30,10,1000,2,10,n,N);
InfoLTplain:=proc(S,K1,K2,n,N) local gu,ope,gadol,K,gu1,i,ope1,lu,j,i1:



ope:=LTrec(convert(S,set),n,N):

gu1:=LTseq(convert(S,set),max(K1,degree(ope,N))):

gadol:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu1)],K2)[K2]:

ope1:=add(coeff(ope,N,i)*mul(n+j,j=1..i)*N^i,i=0..degree(ope,N)):


lu:=As(ope1,n,N):

if lu=FAIL  then
 RETURN([gu1,ope,gadol,FAIL,FAIL ]):
else 
RETURN([gu1,ope,gadol,lu,evalf(lu) ]):
fi:



end:




#InfoLTv(S,K1,K2,n,N): Verbose version of InfoLTv(S,K1,K2,C,n,N) (q.v.) Try:
#InfoLTv({1,2,3},30,10,1000,2,10,n,N);
InfoLTv:=proc(S,K1,K2,n,N) local gu,ope,lu,C1,gadol,i,CONST:
gu:=InfoLT(S,K1,K2,n,N):


print(`On the Number of Labelled Trees where all vertices must have degrees in the set`, {op(S)}):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let a(n) be the the number of labeled trees where all vertices have degrees in the set`, S):
print(``):
print(`Then, for the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(gu[1]):


ope:=gu[2]:

print(`a(n) satsifies the following linear recurrence of order`, degree(ope,N)):
print(``):
print(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0):
print(``):

gadol:=gu[3]:
lu:=gu[4]:

if lu<>FAIL then
print(`a(n) is asymptotically, for some constant, C`):
print( CONST*n!*lu[1]^n*n^lu[2]):
print(``):
print(`and in floating point`):
print(``):
print( CONST*n!*evalf(lu[1])^n*n^lu[2]):
C1:=evalf(gadol/(K2!*evalf(lu[1])^K2*K2^lu[2])):

print(`The constant, CONST, is approximately`, C1):
fi:

print(`Just for fun, the`, K2, `-th term is`):
print(``):
print(gadol):
print(``):
print(`------------------------------------------------------------`):
print(``):
if gu[7]<>FAIL then
print(`Finally, the limiting distribution of the vertices whose degrees are`, S, `are (approximately) respectively`, gu[6]):
print(``):
print(`and the respective limits of the variances (divided by n) are, respectively`, gu[7]):
fi:
end:




#LTpaper(M,K1,n): inputs a positive integer M, larger than 3, a positive integer K1, and a positive integer C and a symbol n gives
#lots of information from InfoLT (q.v.) for all subsets of {1,...,M} with at least three elements and that contain 1. Try:
#LTpaper(3,30,10,n);
LTpaper:=proc(M,K1,n) local i,a,N,gu,gu1,gu2,G,ope,lu,C1,gadol,CONST,L:

if not (type(M,integer) and M>=3) then
print(M, `should be at least 3`):
RETURN(FAIL):
fi:

L:=[]:

print(`--------------------------------------------------------`):
print(`Sequences Enumerating Labeled trees where all vertices are restricted to have degrees in a given set S`):
print(`for all possible subsets of three or more of`, {seq(i,i=1..M)}, `that must include 1 of course`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In the list below we give the first`, K1, `terms , followed by the linear recurrence, followed by the asympotics.`):
print(``):
print(`-------------------------------------------------`):
print(``):
gu:=combinat[powerset]({seq(i,i=2..M)}) minus {{},seq({i},i=2..M)}:

for gu1 in gu do
gu2:=sort(convert(gu1 union {1},list)):
G:=InfoLT(gu2,K1,5000,n,N):
print(`If the set of allowed vertex-degrees is`, gu2, `then the first`, K1, `terms are`):
print(G[1]):

 ope:=G[2]:
print(`The enumerating sequence satisfies the recurrence`):
print(``):
lprint(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0):
print(``):

gadol:=G[3]:

lu:=G[4]:

if lu<>FAIL then
C1:=evalf(gadol/(5000!*evalf(lu[1])^5000*5000^lu[2])):

if (C1>0.01  and C1<100) then

print(`The asymptotics is`):
print(``):

C1:=evalf(gadol/(5000!*evalf(lu[1])^5000*5000^lu[2])):
lprint( CONST*n!*evalf(lu[1])^n*n^lu[2]):
print(``):
print(`where the CONST. is roughly`, evalf(C1,4)):
print(``):

L:=[op(L),[gu1 union {1},evalf(lu[1])]]:
fi:

fi:

print(``):
print(`------------------------------------------------------------`):
print(``):
if G[6]<>FAIL and G[7]<>FAIL  then
print(`Finally, the limiting distribution of the vertices whose degrees are`, gu2, `are (approximately) respectively`, G[6]):
print(``):
print(`and the respective limits of the variances (divided by n) are, respectively`, G[7]):
fi:

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

print(`To sum up, here are all the sets followed by their growth constants`):
print(``):
lprint(L):
print(``):
end:



#LTpaperPlain(M,K1,n): inputs a positive integer M, larger than 3, a positive integer K1, and a positive integer C and a symbol n gives
#lots of information from InfoLTplain (q.v.) for all subsets of {1,...,M} with at least three elements and that contain 1. Try:
#LTpaperPlain(3,30,10,n);
LTpaperPlain:=proc(M,K1,n) local i,a,N,gu,gu1,gu2,G,ope,lu,C1,gadol,CONST,L:

if not (type(M,integer) and M>=3) then
print(M, `should be at least 3`):
RETURN(FAIL):
fi:

L:=[]:

print(`--------------------------------------------------------`):
print(`Sequences Enumerating Labeled trees where all vertices are restricted to have degrees in a given set S`):
print(`for all possible subsets of three or more of`, {seq(i,i=1..M)}, `that must include 1 of course`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In the list below we give the first`, K1, `terms , followed by the linear recurrence, followed by the asympotics.`):
print(``):
print(`-------------------------------------------------`):
print(``):
gu:=combinat[powerset]({seq(i,i=2..M)}) minus {{},seq({i},i=2..M)}:

for gu1 in gu do
gu2:=sort(convert(gu1 union {1},list)):
G:=InfoLTplain(gu2,K1,5000,n,N):
print(`If the set of allowed vertex-degrees is`, gu2, `then the first`, K1, `terms are`):
print(G[1]):

 ope:=G[2]:
print(`The enumerating sequence satisfies the recurrence`):
print(``):
lprint(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0):
print(``):

gadol:=G[3]:

lu:=G[4]:


if lu<>FAIL then
C1:=evalf(gadol/(5000!*evalf(lu[1])^5000*5000^lu[2])):

if (C1>0.01  and C1<100) then

print(`The asymptotics is`):
print(``):

C1:=evalf(gadol/(5000!*evalf(lu[1])^5000*5000^lu[2])):
lprint( CONST*n!*evalf(lu[1])^n*n^lu[2]):
print(``):
print(`where the CONST. is roughly`, evalf(C1,4)):
print(``):

L:=[op(L),[gu1 union {1},evalf(lu[1])]]:
fi:

fi:

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

print(`To sum up, here are all the sets followed by their growth constants`):
print(``):
lprint(L):
print(``):
end:











#LTFpaper(M,k,K,n): inputs a positive integers M, and outputs the first K terms of the sequence enumerating the sequence
#number of labeled trees where no vertex has degree in the set S, for all subsets of {2,3,...,M} with at most k members.
#It also outputs the estimated asympotics of the form n^(n-2)*mu^n. It also outputs the list of the [S,mu] Try:
#LTFpaper(4,1,100,n);
LTFpaper:=proc(M,k,K,n) local S,i1,gu,k1,S1,s1,L,mu,C,gu1,i:
S:={seq(i1,i1=2..M)}:

L:=[]:

print(`Enumerating sequences, and estimated asymptotics, for Enumerating Labeled  trees where no vertex-degrees belong to a prescribed set S, for all subsets of`, {seq(i1,i1=2..M)}, `with at most`, k, ` members `):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

for k1 from 1 to k do
S1:=combinat[choose](S,k1):

for s1 in S1 do
gu:=LTseqF(s1,K):
print(`The first`, K, `terms of the seqeunce enumerating labeled trees where no vertex has a number of neighbors in the set`, s1, ` are `):
print(``):
lprint(gu):
print(``):
gu1:=[seq(gu[i]/i^(i-2),i=1..K)]:
mu:=Zinn(gu1)[2]:
print(``):
C:=evalf(gu[K]/(K^(K-2)*mu^K)):
print(`The asympotics is approximately`):
print(``):
print(C*n^(n-2)*mu^n):
print(``):
L:=[op(L),[s1,mu]]:
od:
od:
print(``):
print(`-----------------------------------------`):
print(``):
print(`To sum up here is a list of the growth constants next to each forbidden set`):
print(``):
lprint(L):

L:
end:

###Start Statistics part




#Start From HISTABRUT

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


#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(FAIL):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

evalf([gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]):

end:



#D1D2(f,x,y,i,j): Given an expression f in x and y, outputs (x*Dx)^i*(y*Dy)^j f . try:
#D1D2(x^4*y^4,x,y,2,3);
D1D2:=proc(f,x,y,i,j) local i1,j1,f1:
f1:=f:
for i1 from 1 to i do
 f1:=expand(x*diff(f1,x)):
od:
for j1 from 1 to j do
 f1:=expand(y*diff(f1,y)):
od:
f1:
end:


#ScaledMixedMoms(f,x,y,N): Given a probability generating function
#f(x,y) (a polynomial, rational function and possibly beyond)
#returns a double L such that
#L[1][1] is the correlation  (in decimals) and L[i][j] is the scaled (i,j) moment
#For example, try:
#ScaledMixedMoments(((x+y+1)^100,x,y,4);

ScaledMixedMoms:=proc(f,x,y,N) local mu,F,memu1,memu2,i,j,T,gu,m1,m2:
mu:=simplify(subs({x=1,y=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,y=1},diff(F,x))):
memu2:=simplify(subs({x=1,y=1},diff(F,y))):

T[1,1]:=1:
T[1,2]:=evalf(memu1):
T[2,1]:=evalf(memu2):

F:=F/(x^memu1*y^memu2):


for i from 3 to N+1 do
T[i,1]:=evalf(subs({x=1,y=1},D1D2(F,x,y,i-1,0))):
od:


for i from 3 to N+1 do
T[1,i]:=evalf(subs({x=1,y=1},D1D2(F,x,y,0,i-1))):
od:


for i from 2 to N+1 do
 for j from 2 to N+1 do
  T[i,j]:=evalf(subs({x=1,y=1},D1D2(F,x,y,i-1,j-1))):
 od:
od:

gu:=[seq([seq(T[i,j],j=1..N+1)],i=1..N+1)]:
m1:=gu[1][2]:
m2:=gu[2][1]:
[seq([seq(gu[i+1][j+1]/(m1^(i/2)*m2^(j/2)),j=1..nops(gu[i+1])-1)],i=1..nops(gu)-1)]:
end:


#Fc(c,x,y): sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2+c*x*y)
Fc:=proc(c,x,y): sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2+c*x*y):end:









#LimNor(c,N): the list of lists of length N, let's call it L, such that its L[i][j] is the scaled momement of the joint normal distribution with
#pdf sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2-c*x*y). Try:
#LimNor(1/2,4);
LimNor:=proc(c,N) local x,y,i,j,f,mu20,mu02:

mu20:=MixMomNorF(c,2,0):
mu02:=MixMomNorF(c,0,2):
[seq([seq(simplify(MixMomNorF(c,i,j)/mu20^(i/2)/mu02^(j/2)),j=1..N)],i=1..N)]:
end:


#LimNor(c,N): the list of lists of length N, let's call it L, such that its L[i][j] is the scaled momement of the joint normal distribution with
#pdf sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2-c*x*y). Try:
#LimNor(1/2,4);
LimNorOld:=proc(c,N) local x,y,i,j,f,mu20,mu02:
f:=Fc(c,x,y):
mu20:=int(int(x^2*f,x=-infinity..infinity),y=-infinity..infinity):
mu02:=int(int(y^2*f,x=-infinity..infinity),y=-infinity..infinity):
[seq([seq(
simplify(int(int(x^i*y^j*f,x=-infinity..infinity),y=-infinity..infinity)/mu20^(i/2)/mu02^(j/2))
,j=1..N)],i=1..N)]:
end:

#End   From HISTABRUT

#LTseqWEslow(P,N,g): The first N terms of the sequence weight-enumerating labeled trees where all the vertex-degrees must belong to the set P according to the weight : product of g[degree(v)] over all vertices v.
#Try:
#LTseqWEslow({1,2,3},30,g);

LTseqWEslow:=proc(P,N,g) local i,L,f:

if not (type(N,posint) and type(P,set) and member(1,P) and {seq(type(i,integer), i in P)}={true}) then
print(`bad input`):
 RETURN(FAIL):
fi:

L:=FEtoSeq(expand(add(g[i]*f^(i-1)/(i-1)!, i in P)),f,N):

L:=expand(SubsFPS(L,add(g[i]*f^i/i!,i in P),f)):
[seq(L[i]*(i-1)!,i=1..N)]:


end:



#LTrecG(P,n,N,g): The linear recurrence operator annihilating the sequence weight-enumerating labelled trees where all vertices have degrees that belong to P, according to the weight Prod(g[degree(v)], over all vertices. Try:
#LTrecG(P,n,N,g);
LTrecG:=proc(P,n,N,g) local z,G,i,ope,c:
G:=add(g[i]*z^(i-1)/(i-1)!, i in P):
ope:=AZd(G^n/z^(n-1)*(n-2)!,z,n,N)[1]:
c:=lcoeff(ope,N):
add(factor(normal(coeff(ope,N,i)/c))*N^i,i=0..degree(ope,N)):
end:



#LTseqWE(P,N,g): A fast version of LTseqWEslow(P,n,g).
#Try:
#LTseqWE({1,2,3},30,g);

LTseqWE:=proc(P,N,g) local ope, INI,n,Sn:
option remember:

if not member(1,P) then
 RETURN(FAIL):
fi:

ope:=LTrecG(P,n,Sn,g):

INI:=LTseqWEslow(P,degree(ope,Sn),g):

expand(SeqFromRec(ope,n,Sn,INI,N)):

end:


#LTstat(P,K1,r):Inputs a LIST,P,  of positive integers (including 1) of length at least 3, and a  large integer K1, and a small integer r, outputs
#the estimated average/n, variance/n, and the scaled 3rd through r-th moment of the random variables "number of vertices of degree P[i]"
#among labeled trees with n vertices, where all the vertices must be of degree in P
#where the average and variance are divided by n. For n from K1-20 to K1. Try:
#LTstat([1,2,3],1000,6);
LTstat:=proc(P,K1,r) local i,L,g,L1,M,M1,M11,i1,err1,err2,d1,d2:

if not (type(P,list) and {seq(type(P[i],posint),i=1..nops(P))}={true} and type(r,posint) and r>=2 and type(K1,posint) and K1>=30 and nops(P)>=3 ) then
 print(`bad input`):
 RETURN(FAIL):
fi:

L:=LTseqWE(convert(P,set),K1,g):

M:=[]:

for i from 1 to nops(P) do
 L1:=subs({seq(g[P[i1]]=1,i1=1..i-1),seq(g[P[i1]]=1,i1=i+1..nops(P))},L):
  M1:=[]:
   for i1 from K1-1 to K1 do
    M11:=Alpha(L1[i1],g[P[i]],r):
    if M11=FAIL then
       RETURN(FAIL):
    fi:
    M11:=[M11[1]/i1,M11[2]/i1,op(3..r,M11)]:
    M1:=[op(M1),M11]:
  od:

err1:=M1[1][1]-M1[2][1]:
err2:=M1[1][2]-M1[2][2]:
d1:=trunc(-log[10](abs(err1))):
d2:=trunc(-log[10](abs(err2))):

M:=[op(M),[evalf(M1[2][1],d1), evalf(M1[2][2],d2), seq(M1[2][i1],i1=3..r)] ]:
od:
M:
end:




#LTstatJ1(P,K1,r,deg1,deg2):Inputs a LIST of positive integers (including 1) and a  large integer K1, and a small integer r, 
#and two members deg1,deg2, of P outputs scaled mixed moments for the random variables "number of vertices with degree deg1
#and "number of vertices of degree deg2" up to the r-th mixed moment. Try:
#LTstatJ1([1,2,3],1000,2,1,2);
LTstatJ1:=proc(P,K1,r,deg1, deg2) local i,L,g,L1,M1,M11,i1,err1,d1:

if not (type(P,list) and {seq(type(P[i],posint),i=1..nops(P))}={true} and type(r,posint) and r>=2 and type(K1,posint) and K1>=30 and member(deg1,{op(P)}) and member(deg2, {op(P)}) and deg1<deg2 ) then
 print(`bad input`):
 RETURN(FAIL):
fi:

L:=LTseqWE(convert(P,set),K1,g):


 L1:=subs({seq(g[P[i1]]=1,i1=1..deg1-1),seq(g[P[i1]]=1,i1=deg1+1..deg2-1), seq(g[P[i1]]=1,i1=deg2+1..nops(P))},L):


  M1:=[]:
   for i1 from K1-1 to K1 do
    M11:=ScaledMixedMoms(L1[i1],g[deg1],g[deg2],r):
       M1:=[op(M1),M11]:
  od:


err1:=M1[1][1][1]-M1[2][1][1]:

d1:=trunc(-log[10](abs(err1))):
evalf(M1[2],d1):


end:




#LTseqFwe(P,N,g): The first N terms of the sequence weight-enumerating labeled trees according to the weight Prod(g[i]^(Number of vertices with degree i) in the list of psoitive integers P
#Try:
#LTseqFwe([1,2],10,g);
LTseqFwe:=proc(P,N,g) local i,L,f:

if not (type(N,posint) and type(P,list) and {seq(type(i,integer), i in {op(P)})}={true}) then
print(`bad input`):
 RETURN(FAIL):
fi:


L:=FEtoSeq(exp(f)+add((g[i]-1)*f^(i-1)/(i-1)!, i in P),f,N):
L:=SubsFPS(L,exp(f)+add((g[i]-1)*f^i/i!,i in P),f):
[seq(L[i]*(i-1)!,i=1..N)]:


end:


#LTseqFwe1(P,n,g): The  weight-enumerator  labeled trees  with exactlty n vertices using the Lagrange inversion formula according to the weight Prod(g[i]^(Number of vertices with degree i) in the list pf psoitive integers P
#Try:
#LTseqFwe1([1,2],10,g);
LTseqFwe1:=proc(P,n,g) local i,f,z:
if n=1 then
 RETURN(1):
fi:

f:=exp(z)+add((g[P[i]]-1)*z^(P[i]-1)/(P[i]-1)! ,i=1..nops(P)):

expand((n-2)!*coeff(taylor(f^n,z=0,n+4),z,n-2)):

end:



#LTmom(P,n,a): inputs a list of positive integers P, a symbol n, and a list of positive integers a (of the same length of P) outputs the
#mixed-a-moments of the random variables (in the population of labelled trees with n vertices) "number of ocurrences of degree P[i] vertices/
#Try:
#LTmom([1],n,[2]);
LTmom:=proc(P,n,a) local i,j,f,z,g,Z,lu,lu1,lu11:
option rememer:
if not (type(P,list) and type(a,list) and nops(P)=nops(a) and type(n,symbol)) then
print(`Bad input`):
RETURN(FAIL):
fi:


f:=exp(z)+add((g[P[i]]-1)*z^(P[i]-1)/(P[i]-1)! ,i=1..nops(P)):

f:=f^n:

for i from 1 to nops(P) do
for j from 1 to a[i] do
 f:=simplify(g[P[i]]*diff(f,g[P[i]])):
od:
od:

for i from 1 to nops(P) do
f:=subs(g[P[i]]=1,f):
od:

f:=expand(f):
f:=subs(exp(z)=Z,f):

f:=simplify(f/Z^n):
f:

f:=expand(f):


lu:=0:

for i from 0 to -ldegree(f,Z) do
 lu1:=coeff(f,Z,-i):
  for j from 0 to degree(lu1,z) do
    lu11:=coeff(lu1,z,j):

 lu:=lu+lu11*(n-i)^(n-2-j)*simplify((n-2)!/(n-2-j)!):
od:
od:

lu/n^(n-2):

end:

#LTmom1(d,n,a): The a-th moment (not about the mean) of the random variable "number of vertices of degree d" in the set of labeled trees on n vertices. Try:
#LTmom1(1,n,2);
LTmom1:=proc(d,n,a)
option remember:
LTmom([d],n,[a]):
end:



#LTmom1am(d,n,a): The a-th moment (about the mean) of the random variable "number of vertices of degree d" in the set of labeled trees on n vertices. Try:
#LTmom1am(1,n,2);
LTmom1am:=proc(d,n,a)  local mu,k:
mu:=LTmom1(d,n,1):
expand(add(binomial(a,k)*(-mu)^(a-k)*LTmom1(d,n,k),k=0..a)  ):
end:



#LTmom1amL(d,a): The limit of the a-th moment (about the mean) of the random variable "number of vertices of degree d" divided by n^(trunc(d/2)) in the set of labeled trees on n vertices. Try:
#LTmom1amL(1,n,2);
LTmom1amL:=proc(d,a)  local n,gu:
gu:=LTmom1am(d,n,a):
limit(gu/n^(trunc(a/2)),n=infinity):
end:






#LTmom2(d1,d2,n,a1,a2): The (a1,a2) mixed moments (not about the mean) of the random variables "number of vertices of degree d1"  abd "number of vertices with degree d2", in the set of labeled trees on n vertices. Try:
#LTmom2(1,2,n,2,2);
LTmom2:=proc(d1,d2,n,a1,a2)
option remember:
LTmom([d1,d2],n,[a1,a2]):
end:


#LTmom2am(d1,d2,n,a1,a2): The (a1,a2) mixed moment (about the mean) of the random variables "number of vertices of degree d1"  and "number of vertices of degree d2" in the set of labeled trees on n vertices. Try:
#LTmom2am(1,2,n,2,2);
LTmom2am:=proc(d1,d2,n,a1,a2)  local mu1,mu2, k1,k2:
mu1:=LTmom1(d1,n,1):
mu2:=LTmom1(d2,n,1):
expand(add(add(binomial(a1,k1)*(-mu1)^(a1-k1)*binomial(a2,k2)*(-mu2)^(a2-k2)*LTmom2(d1,d2,n,k1,k2),k1=0..a1),k2=0..a2)  ):
end:


#LTmom2amL(d1,d2,n,a1,a2): The Limit of LTtom2(d1,d2.n,a1,a2)  over (n^(a1+a2)/2). Try:
#LTmom2amL(1,2,n,2,2);
LTmom2amL:=proc(d1,d2,a1,a2)  local n:
limit(LTmom2am(d1,d2,n,a1,a2)/n^((a1+a2)/2),n=infinity):
end:



#LimMixedMoms(d1,d2,L): The limit of the scaled momemts of the random variables "number of vertices with degree d1" and "number of vertices of degree d2" in labeled trees with n vertices, as n goes to infinity. Try
#LimMixedMoms(1,2,2); 
LimMixedMoms:=proc(d1,d2,L) local j1,j2,V1,V2,n,gu1:
option remember:
V1:=limit(LTmom1am(d1,n,2)/n,n=infinity):
V2:=limit(LTmom1am(d2,n,2)/n,n=infinity):
gu1:=simplify([seq([seq(LTmom2amL(d1,d2,j1,j2)/(V1^(j1/2)*V2^(j2/2)),j2=1..L)],j1=1..L)]):
gu1:
end;



#AveFor(d,n): The rormula for the average number of degree-d vertices in a labeled tree on n vertices. Try:
#AbeFor(d,n);
AveFor:=proc(d,n):simplify(n!/(n-d-1)!*(n-1)^(n-d-2)/(d-1)!/n^(n-2)):end:

#LimAveFor(d): the formula ,exp(-1)/(d-1)!, for the limit of the expectation of the r.v. "Number of verticecs in a random labeled tree on n vertices with degree d" divided by n. Try:
#LimAveFor(d);
LimAveFor:=proc(d): exp(-1)/(d-1)!:end:


#LimVarFor(d): the formula ,exp(-1)/(d-1)!- (d^2-4*d+5)/(d-1)!^2, for the limit of the variance of the r.v. "Number of verticecs in a random labeled tree on n vertices with degree d" divided by n. Try:
#LimVarFor(d);
LimVarFor:=proc(d): 
exp(-1)/(d-1)!- exp(-2)*(d^2-4*d+5)/(d-1)!^2:
end:



#LimCovFor(d1,d2): the formula ,exp(-2)*(2*d1+2*d2-d1*d2-5)/((d1-1)!*(d2-1)!), for the limit of the covariance of the r.v.s "Number of verticecs in a random labeled tree on n vertices with degrees d1 and d2" divided by n. Try:
#LimCovFor(d1,d2);
LimCovFor:=proc(d1,d2): 
exp(-2)*(2*d1+2*d2-d1*d2-5)/((d1-1)!*(d2-1)!):
end:

#LimCorFor(d1,d2): the formula for the limit of the correlation of the r.v.s "Number of verticecs in a random labeled tree on n vertices with degrees d1 and d2" divided by n. Try:
#LimCorFor(d1,d2);
LimCorFor:=proc(d1,d2):
LimCovFor(d1,d2)/sqrt(LimVarFor(d1)*LimVarFor(d2)):
end:


##start checking procedures

#DegSeq(G,n) : Given a graph G on n vertices, outputs the list of length n such that its i-th entry is the number of vertices with degree i. Try
#DegSeq(GG(5),5);
DegSeq:=proc(G,n) local i,T,e:
for i from 1 to n do
 T[i]:=0:
od:

for e in G do
 T[e[1]]:=T[e[1]]+1:
 T[e[2]]:=T[e[2]]+1:
od:

[seq(T[i],i=1..n)]:
end: 


#Wt(G,g): The weight of the graph G according to product(vertices) g[i]^degree[i]
Wt:=proc(G,n,g) local i,L:
L:=DegSeq(G,n):
mul(g[L[i]],i=1..n):
end:


#CheckLTseqWE(P,N): Checks LTseqWE(P,N,g).
#Try:
#CheckLTseqWE({1,2,3},8);
CheckLTseqWE:=proc(P,N) local g,gu,n1,L,L1,gu1:
L:=LTseqWE(P,N,g) :

L1:=[]:

for n1 from 1 to N do
 gu:=STp(n1,CG(n1),P):
 L1:=[op(L1),add(Wt(gu1,n1,g),gu1 in gu)]:
od:

evalb(L=L1):

end:



#CheckLTseqFwe(P,N): Checks LTseqFwe(P,N,g).
#Try:
#CheckLTseqFwe({1,2,3},8);
CheckLTseqFwe:=proc(P,N) local g,gu,n1,L,L1,gu1,j,lu:
L:=LTseqFwe(P,N,g):
L1:=[1]:

for n1 from 2 to N do
 gu:=ST(n1,CG(n1)):
 lu:=add(Wt(gu1,n1,g),gu1 in gu):
 for j from 1 to N do
  if not member(j,{op(P)}) then
   lu:=subs(g[j]=1,lu):
  fi:
od:
 L1:=[op(L1),lu]:

od:

evalb(L=L1):

end:





#CheckLTseq: Checks empirically LTseq. Try:
#CheckLTseq({1,2,3},6);
CheckLTseq:=proc(P,N)  local n1:
evalb(LTseq(P,N)=[seq(nops(STp(n1,CG(n1),P)),n1=1..N)]):
end:




#CheckLTseq: Checks empirically LTseq. Try:
#CheckLTseq({1,2,3},6);
CheckLTseq:=proc(P,N)  local n1:
evalb(LTseq(P,N)=[seq(nops(STp(n1,CG(n1),P)),n1=1..N)]):
end:


#CheckLTseqF: Checks empirically LTseqF. Try:
#CheckLTseqF({2,3},6);
CheckLTseqF:=proc(P,N) local n1:
if member(1,P) then
 print(`1 can't be a member`):
 RETURN(FAIL):
fi:
evalb(LTseqF(P,N)=[seq(nops(STf(n1,CG(n1),P)),n1=1..N)]):
end:



#CheckLTmom1(d,N,a): checks LTmom1(d,n,a) for n=1..N. Try:
#CheckLTmom1(1,6,1);
CheckLTmom1:=proc(d,N,a) local n1,g,mu,lu1,n,gu,lu:

gu:=LTmom1(d,n,a):

for n1 from a+1 to N do
 
 lu:=ST(n1,CG(n1)):

mu:=0: 
for lu1 in lu do
  mu:=mu+degree(Wt(lu1,n1,g),g[d])^a:
 od:
 mu:=mu/nops(lu):
 if mu<>subs(n=n1,gu) then
   RETURN(false):
 fi:
od:
true:
end:




#CheckLTmom1am(d,N,a): checks LTmom1am(d,n,a) for n=1..N. Try:
#CheckLTmom1am(1,6,1);
CheckLTmom1am:=proc(d,N,a) local n1,g,mu,lu1,n,gu,lu,ave:

gu:=LTmom1am(d,n,a):
ave:=LTmom1(d,n,1):
for n1 from a+1 to N do
 
 lu:=ST(n1,CG(n1)):

mu:=0: 
for lu1 in lu do
  mu:=mu+(degree(Wt(lu1,n1,g),g[d])-subs(n=n1,ave))^a:
 od:
 mu:=mu/nops(lu):
 if mu<>subs(n=n1,gu) then
   RETURN(false):
 fi:
od:
true:
end:





##End checking procedures

#T1rn(r,n): the number of labeled trees on n vertices where none of the vertices have degree r. Try:
#T1rn(2,20);
T1rn:=proc(r,n) local k,s:

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

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

s:=r-1:
add((-1)^k/s!^k*binomial(n,k)*(n-k)^(n-2-k*s)*(n-2)!/(n-2-k*s)!,k=0..trunc((n-2)/s)):
end:



#T2rn(r1,r2,n): the number of labeled trees on n vertices where none of the vertices  have degrees that belong to {r1,r2}. Try:
#T2rn(2,3,20);
T2rn:=proc(r1,r2,n) local k1,k2,s1,s2:
if not (type(r1,integer) and type(r2,integer) and min(r1,r2)>1) then
 RETURN(FAIL):
fi:
if n=1 then
 RETURN(0):
fi:
s1:=r1-1:
s2:=r2-1:
add(
add((-1)^(k1+k2)/(s1!^k1*s2!^k2)*n!/(k1!*k2!*(n-k1-k2)!)*(n-k1-k2)^(n-2-k1*s1-k2*s2)*(n-2)!/(n-2-k1*s1-k2*s2)!,k2=0..trunc((n-2-s1*k1)/s2)),
k1=0..trunc((n-2)/s1)):
end:



#T3rn(r1,r2,r3,n): the number of labeled trees on n vertices where none of the vertices  have degrees that belong to {r1,r2,r3}. Try:
#T3rn(2,3,4,20);
T3rn:=proc(r1,r2,r3,n) local k1,k2,k3,s1,s2,s3:
if not (type(r1,integer) and type(r2,integer) and type(r3,integer) and min(r1,r2,r3)>1) then
 RETURN(FAIL):
fi:
if n=1 then
 RETURN(0):
fi:
s1:=r1-1:
s2:=r2-1:
s3:=r3-1:

add(
 add(
  add(
     (-1)^(k1+k2+k3)/(s1!^k1*s2!^k2*s3!^k3)*n!/(k1!*k2!*k3!*(n-k1-k2-k3)!)*(n-k1-k2-k3)^(n-2-k1*s1-k2*s2-k3*s3)*(n-2)!/(n-2-k1*s1-k2*s2-k3*s3)!,
       k3=0..trunc(n-2-s1*k1-s2*k2)/s3),
    k2=0..trunc((n-2-s1*k1)/s2)    ),
k1=0..trunc((n-2)/s1
)
):

end:






#start using recurrences for the mixed moments of the bivariate normal distribution
#Mixed moments  of the bivariate distribution with collrelation (2*m,2*n)
NorE:=proc(c,m,n) option remember:

if m=0  then

if n=0 then
 RETURN(1):
 
elif n=1 then

RETURN(1/(1-c^2)):

else

RETURN(normal(-(2*c^2*m-2*c^2*n+2*c^2+4*n-3)/(c-1)/(c+1)*NorE(c,m,n-1)+2*(n-1)*(2*n-3)/(c-1)/(c+1)*NorE(c,m,n-2))):
fi:

elif m=1  then

if n=0 then

RETURN(1/(1-c^2)):

elif n=1 then

RETURN((1+2*c^2)/(1-c^2)^2):

else

RETURN(normal(-(2*c^2*m-2*c^2*n+2*c^2+4*n-3)/(c-1)/(c+1)*NorE(c,m,n-1)+2*(n-1)*(2*n-3)/(c-1)/(c+1)*NorE(c,m,n-2))):
fi:

else

RETURN( normal((2*c^2*m-2*c^2*n-2*c^2-4*m+3)/(c-1)/(c+1)*NorE(c,m-1,n)+2*(m-1)*(2*m-3)/(c-1)/(c+1)*NorE(c,m-2,n))):
fi:
end:

#Mixed moments  of the bivariate distribution with collrelation (2*m+1,2*n+1)
NorO:=proc(c,m,n) option remember:

if m=0  then

if n=0 then
 RETURN(c/(1-c^2)):
 
elif n=1 then

RETURN(3*c/(1-c^2)^2):

else


RETURN(normal(-(2*c^2*m-2*c^2*n+2*c^2+4*n-1)/(c-1)/(c+1)*NorO(c,m,n-1)+2*(2*n-1)*(n-1)/(c-1)/(c+1)*NorO(c,m,n-2))):
fi:

elif m=1  then

if n=0 then

RETURN(3*c/(1-c^2)^2):

elif n=1 then

RETURN(3*c*(3+2*c^2)/(1-c^2)^3):

else
RETURN(normal(-(2*c^2*m-2*c^2*n+2*c^2+4*n-1)/(c-1)/(c+1)*NorO(c,m,n-1)+2*(2*n-1)*(n-1)/(c-1)/(c+1)*NorO(c,m,n-2))):

fi:

else

RETURN(normal((2*c^2*m-2*c^2*n-2*c^2-4*m+1)/(c-1)/(c+1)*NorO(c,m-1,n)+2*(2*m-1)*(m-1)/(c-1)/(c+1)*NorO(c,m-2,n))):
fi:
end:


#MixMomNorF(c,m,n): The mixed-moment of the bivairate normal distribution with correlation c
#Try: MixNorF(c,4,6);
MixMomNorF:=proc(c,m,n)
if m+n mod 2=1 then
 RETURN(0):
elif m mod 2=0 and n mod 2=0 then
 RETURN(NorE(c,m/2,n/2)):
else
 RETURN(NorO(c,(m-1)/2,(n-1)/2)):
fi:
end:

##end using recurrences for the mixed moments of the bivariate normal distribution



#ProveAsyNor(d1,r): Rigorously proves the asymptotic normality of the random varialbe "number of vertices with degree d" using the method of moments by checking that the scaled ven moments are the same as those of the normal distribution
#i.e. 1,3,15,105, ... to the 2r-th moment. Try:
#ProveAsyNor(1,5);
ProveAsyNor:=proc(d1,r) local i,m2:

m2:=LTmom1amL(d1,2):
evalb([seq(simplify(LTmom1amL(d1,2*i)/m2^i),i=1..r)]=[seq((2*i)!/(i!*2^i),i=1..r)]):

end:

#ProveJointAsyNor(d1d2,L): Rigorously proves the joint asymptotic (but, not independent) normality of the random varialbes 
#"number of vertices with degree d1"  and "number of vertices with degree d2" in labeled trees
#using the method of moments by checking that the scaled mixed moments are the same as those of the scaled mixed 
#moments with the same correlation. Try
#ProveJointAsyNor(1,2,4);
ProveJointAsyNor:=proc(d1,d2,L)  local gu1,gu2,j1,j2:
gu1:=LimMixedMoms(d1,d2,L):
gu2:=simplify(LimNor(gu1[1][1],L)):
evalb({seq(seq(simplify(gu1[j1][j2]-gu2[j1][j2]),j1=1..L),j2=1..L)}={0}):
end:



#StatPaper(MaxDeg,MaxMom1, MaxMom2): Inputs positive integers MaxDeg, MaxMom1, and MaxMom2 and outputs a paper about the r.v. "number of vertices of degree d"
#with explicit expressions for the moments for d between 1 amd MaxMom1 and for the mixed moments up to  MaxMom2. Try:
#StatPaper(4,4,4);
StatPaper:=proc(MaxDeg,MaxMom1, MaxMom2) local n,X,d,lu,d1,d2,v2,m1,m2,i1:
print(``):
print(`Statistical Analysis of the Random Variables "Number of Vertices with d Neighbors", in the Set of labelled trees,for genreal d but especially for d between 1 and`, MaxDeg, `and Explicit expressions for their moments up to`, MaxMom1):
print(`and mixed moments up to `, MaxMom2):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In this article we will study the random variable let's call it`, X[d], `that is the number of vertices in  a random labelled tree on n vertices that have d neighbors. In particular, X[1], is the number of leaves`):
print(`Let's start with some general formulas valid for all degrees `):
print(``):
print(`Theorem 1: For all positive integers d the asymptotic expression for the expectation of X[d] is`):
print(``):
lu:=LimAveFor(d):
print(``):
print(lu*n):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu*n):
print(``):
lu:=LimVarFor(d):

print(`The variance of X[d] is, asymptotically`):
print(``):
print(lu*n):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu*n):
print(``):


lu:=LimCovFor(d1,d2):

print(`The covariance of X[d1] and X[d2] is, asymptotically`):
print(``):
print(lu*n):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu*n):
print(``):


lu:=LimCorFor(d1,d2):

print(`Hence The limit of the correlation of X[d1] and X[d2], as n goes to infinity  is`):
print(``):
print(lu):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):
print(`Now let's give explicit formulas for all the moments, about the mean, of X[d] from d=1 to`, d=MaxDeg):
print(`For the second throgh the`, MaxMom1, ` moment `):

print(``):

for d1 from  1 to MaxDeg do
print(``):
print(`-------------------------------------------------------------------------------`):
print(``):
for m1 from 2 to MaxMom1 do
print(``):
if m1=2 then
print(`The variance of`, X[d1] , `is `):
else
print(` The `, m1, `-th moment about the mean of`, X[d1] , `is `):
fi:
lu:=LTmom1am(d1,n,m1):
print(``):
print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):
od:

print(`Also the limits of the scaled moments about the mean up to the`, 2*MaxMom1, `are `):
v2:=LTmom1amL(d1,2):
print(``):
print([seq(simplify(LTmom1amL(d1,2*i1)/v2^i1),i1=1..MaxMom1)]):
print(``):
print(`that are the same as those of the Normal distribution, proving asympotic normality up to the`, 2*MaxMom1, `moment. `):
od:
print(``):
print(`-----------------------------------------------------------`):
print(``):

for d1 from 1 to MaxDeg do
 for d2 from d1+1 to MaxDeg do
 print(`We will now study how `, X[d1], `and `, X[d2], `intereact `):

 for m1 from 1 to MaxMom2 do
 for m2 from 1 to MaxMom2 do
  
 lu:=LTmom2am(d1,d2,n,m1,m2):

print(` The `, (m1,m2), `mixed moment of`, X[d1], `and `, X[d2], `is `):
print(``):
print(lu):
print(``):
print(`and in Maple notation `):
print(``):
lprint(lu):
print(``):
od:
od:
print(`In particular we have an indication that`, X[d1], `and `, X[d2], `are joint-asymptoically normal with limiting correlation`):
print(``):
print(LimCorFor(d1,d2)):
print(``):
print(`and in decomals`):

print(``):
lprint(evalf(LimCorFor(d1,d2))):
print(``):
print(`------------------------------------------------`):
od:
od:
print(``):
print(`---------------------------------------------------`):
print(``):
print(`This ends this paper that took`, time(), `seconds to produce. `):

end: