######################################################################
##RPS: Save this file as  RPS                                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read RPS<Enter>                                    #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec. 9, 2011
 
print(`Created:  Dec. 9, 2011`):
print(` This is RPS `):
print(`Named after Richard Peter Stanley`):
print(`It is one of the packages that accompany the article `):
print(`Automatic Solution of Richard Stanley's Amer. Math. Monthly Problem #11610`):
print(`and ANY Problem of That Type`):

print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`available from http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/rps.html`):

print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra11:=proc()

if args=NULL then
 print(` The secondarysupporting procedures are: HSP `):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: algtodiff, algtodiff1, AsyC `):
 print(`AZc, CheckAsy, CheckGWtE, CheckRec,`):
 print(` DiffToRec1, Empir, EqsVars,  findrec, Findrec,`):
 print(` Images1, Khalek, Milim, SeqFromRec, Wt, WtEd, Zugot`):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` ALGg, DOaz, GWtE, Info1, Info2, Mishpat`):
 print(` ProveGF, RECaz, RECg, RSct, RSseq, TerseBook, VerboseBook`):
 print(` WtE `):


elif nops([args])=1 and op(1,[args])=ALGg then
print(`ALGg(m,w1,w2,K,x,P): Given a positive integer`):
print(`m and two words (lists), w1, w2 of the same length in the`):
print(`alphabet {1, ...,m}, `):
print(`and a positive integer K, and symbols x, P`):
print(`uses guessing with K data points`):
print(`to conjecture a polynomial in x and P`):
print(`such that you get zero when`):
print(`you plug-in the ordinary generating function of`):
print(`  the sequence`):
print(`"number of words in {1, ...m} of length n that has`):
print(`as many factors (substrings) that are w1 and`):
print(`factors that are w2 equals `):
print(`For example, to get a recurrence operator for`):
print(`Richard Stanley's`):
print(`Amer. Math. Monthly problem #11610 (page 937`):
dprint(` of v. 118 (No. 10), [Dec. 2011]), please type:`):
print(`ALGg(2,[1,1],[1,2],40,x,P);`):

elif nops([args])=1 and op(1,[args])=algtodiff then
print(`algtodiff(F,P,x,D1): Inputs a polynomial F,`):
print(`in P and x, and outputs`):
print(`a differential operator annihilating`):
print(`P(x) where P(x) is the alg. formal power series`):
print(`satisfying the algebraic equation`):
print(`F(P(x),x))=0 . For example, try:`):
print(`algtodiff(P-1-x*P^2,P,x,D1);`):

elif nops([args])=1 and op(1,[args])=algtodiff1 then
print(`algtodiff1(F,P,x,D1,ORDER): Inputs a polynomial F,`):
print(`in P and x, and an ORDER, and outputs`):
print(`a differential operator annihilating`):
print(`P(x) where P(x) is the alg. formal power series`):
print(`satisfying the algebraic equation`):
print(`F(P(x),x))=0 . For example, try:`):
print(`algtodiff1(P-1-x*P^2,P,x,D1,2);`):


elif nops([args])=1 and op(1,[args])=AsyC then
print(`AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0, with the given initial conditions`):
print(`where ope(n,N) is  a recurrence operator`):
print(`up to the M's term,`):
print(`and complete with an empirically derived constant in front`):
print(`using K terms`):

elif nops([args])=1 and op(1,[args])=AZc then
print(` AZc(A,y,x,D): The Almkvist-Zeilberger Continuous algorithm`):
print(`Given an hyperexponential expression A, of y and x`):
print(`returns a differential operator ope(x,D) that annihilates`):
print(` the function x->Integral( A(x,y),y);`):
print(`For example, try:`):
print(`AZc(1/(1-x*(y+1/y))/y,y,x,Dx);`):

elif nops([args])=1 and op(1,[args])=CheckAsy then
print(`CheckAsy(a,n,L,eps): checks whether the asympotic formula  a, an`):
print(`expression in n, is a good one for the sequence L. With tolerance`):
print(`eps. For example, try:`):
print(`CheckAsy(2^n,n,[seq(2^i,i=1..100)],1/100);`):

elif nops([args])=1 and op(1,[args])=CheckGWtE then
print(`CheckGWtE(m,k,N): checks the first N terms of`):
print(`the output of CheckGWtE(m,k,z,t); For example, try:`):
print(`CheckGWtE(2,2,5);`):

elif nops([args])=1 and op(1,[args])=CheckRec then
print(`CheckRec(L,ope,n,N): Given a sequence of numbers L and a`):
print(`recurrence operator ope(n,N) checks whether indeed`):
print(`L is annihilated by ope(n,N). For example, try:`):
print(`CheckRec([1,2,4,8,16],N-2,n,N);`):

elif nops([args])=1 and op(1,[args])=DiffToRec1 then
print(`DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator`):
print(`ope(t,Dt) and symbols n,N, outputs the linear recurrence`):
print(`operator with polynomial coefficients satisfied by`):
print(`the sequence of coeff. of any f.p.s. annihilating ope(t,Dt)`):
print(`where N is the shift in n. For example, try:`):
print(`DiffToRec1(Dt+t, t,Dt,n,N);`):

elif nops([args])=1 and op(1,[args])=DOaz then
print(`DOaz(m,w1,w2,z,Dz): Given a positive integer`):
print(`m and two words (lists) of the same length, w1, w2 in the`):
print(`alphabet {1, ...,m}, `):
print(`uses the continuous Almkvist-Zeilberger algorithm,`):
print(`to find a linear differential operator with`):
print(`polynomial coefficients, in terms of z and the`):
print(`differentiation w.r.t. z, Dz, that annihilates`):
print(`the (ordinary) generating function`):
print(`of the sequence`):
print(`"number of words in {1, ...m} of length n that has`):
print(`as many factors (substrings) that are w1 as it has`):
print(`factors that are w2. `):
print(``):
print(`For example, to get the differential operator for`):
print(`the generating function for the sequence described in`):
print(`Richard Stanley's `):
print(`Amer. Math. Monthly problem #11610 (page 937`):
print(` of v. 118 (No. 10), [Dec. 2011]), please type:`):
print(`DOaz(2,[1,1],[1,2],t,Dt);`):


elif nops([args])=1 and op(1,[args])=Empir then
print(`Empir(L,x,P): inputs a sequence L and `):
print(`and conjectures an polynomial equation`):
print(`F(x,P)=0 satisfies by the formal power`):
print(`series L[1]+L[2]*x+L[3]*x^2+...`):
print(`For example, try:`):
print(`Empir([seq(i1,i1=1..40)],x,P);`):

elif nops([args])=1 and op(1,[args])=EqsVars then
print(`EqsVars(m,k,z,t,A):The set of equations, followed by the`):
print(`set of variables that enables to compute the generating`):
print(`function in z and t[w1] (w1 in [1,...,m]^k)`):
print(`that is the weight-enumerator of all words in`):
print(`the alphabet {1, ...,m} with weight z^length(w) times`):
print(`product of t[w_i,...,w_i+k-1] for i from 1 to nops(w)-k+1`):
print(`A is a symbol for the indexing of the weight enumerator`):
print(`of words,`):
print(`For example, try`):
print(`EqsVars(2,2,z,t,A);`):

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 nops([args])=1 and op(1,[args])=HSP then
print(`HSP(ope,n,N): the highest singularity of the recurrence`):
print(`operator ope(n,N). For example, try:`):
print(`HSP((n-3)*(n-4)*N+1,n,N);`):

elif nops([args])=1 and op(1,[args])=Images1 then
print(`Images1(m,S): Given a set of words S in the alphabet`):
print(`{1, ...,m}, finds the set of images under permuting`):
print(`the letters and reversing the words.`):
print(`For example, try:`):
print(`Images1(2,{[1,1],[1,2]});`):


elif nops([args])=1 and op(1,[args])=Info1 then
print(`Info1(m,w1,w2,n,N,P,t,K): Given a pos. integer m,`):
print(`and two words w1 and w2 of the same length in the`):
print(`alphabet {1, ..., m}, and symbols n and N`):
print(`P and t.  It also inputs a positive integer`):
print(`K (the number of terms for guessing the alg. relation, see below)`):
print(`We are interested in the integer sequence`):
print(`a(n):=number of n-lettered words in the alphabet {1, ...,m}`):
print(`that has as many substrings (alias factors, `):
print(`alias consecutive subwords) that are w1 as there are w2.`):
print(`The output this procedure  is a list consisting`):
print(`(i) a (conjectured) recurrence operator ope(n,N) annihilating`):
print(`the sequence a(n) together with the initial values `):
print(`(starting at n=1) providing`):
print(`a full holonomic (P-recursive) description of the sequence`):
print(`(ii) The minimal polynomial F(P,x) such that`):
print(`F(P(x),x)=0 where P(x) is the ordinary generating function of`):
print(`the sequence: i.e.,`):
print(`P(x)=Sum(a[n]*x^n,n=0..infinity);`):
print(`For example, for Richard Stanley's Amer. Math. Monthly Problem`):
print(`11610 (118(10)[Dec. 2011], p. 937)`):
print(`type`):
print(`Info1(2,[1,1],[1,2],n,N,P,t,40);`):


elif nops([args])=1 and op(1,[args])=Info2 then
print(`Info2(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m,`):
print(`and two words w1 and w2 of the same length in the`):
print(`alphabet {1, ..., m}, and symbols n and N`):
print(`P and t.  It also inputs a positive integer`):
print(`K (the number of terms for guessing the alg. relation, see below)`):
print(`We are interested in the integer sequence`):
print(`a(n):=number of n-lettered words in the alphabet {1, ...,m}`):
print(`that has as many substrings (alias factors, `):
print(`alias consecutive subwords) that are w1 as there are w2.`):
print(`The output this procedure  is a list consisting`):
print(`(i) a (conjectured) recurrence operator ope(n,N) annihilating`):
print(`the sequence a(n) together with the initial values `):
print(`(starting at n=1) providing`):
print(`a full holonomic (P-recursive) description of the sequence`):
print(`(ii) The minimal polynomial F(P,x) such that`):
print(`F(P(x),x)=0 where P(x) is the ordinary generating function of`):
print(`the sequence: i.e.,`):
print(`P(x)=Sum(a[n]*x^n,n=0..infinity);`):
print(`(iii) an asymptotic expression for a(n) to order M`):
print(`(iv): the ratios of the sequence to the asymptotic expression`):
print(`for the terms 1990 to 2000 (just for checking)`):
print(`The sequene of ratios of the members of the sequence to the asymptotics`):
print(`only for checking purposes`):
print(`(vi): a conjectured value for the constant, if avaiable, or else 0`):
print(`For example, for Richard Stanley's Amer. Math. Monthly Problem`):
print(`11610 (118(10)[Dec. 2011], p. 937)`):
print(`type`):
print(`Info2(2,[1,1],[1,2],n,N,P,t,40,3);`):

elif nops([args])=1 and op(1,[args])=Khalek then
print(`Khalek(ope1,ope2,n,N): given linear recurrence operators ope1`):
print(`and (monin in N) ope2 in n and N (the forward shift operator in N),`):
print(`outputs operators Q(N,n) and R(N,n) such that`):
print(`ope1(N)=Q(N,n)ope2(n,N)+R(N,n) `):
print(` where the order of Q(N,n)(i.e. degree in N)`):
print(`is the difference in the oders of ope1 and ope2 and the order of`):
print(`R(N,n) is less than the order of ope2.`):
print(`For example, try: `):
print(` Khalek(n*N^2+n/(n+1)*N,N+1/n,n,N); `):

elif nops([args])=1 and op(1,[args])=GWtE then
print(`GWtE(m,k,z,t): The global weight-enumerator of words in`):
print(`the alphabet {1,..,m} according to the weight z^length`):
print(`times the product of t[FactorsOfLength_k]. For`):
print(`example, try:`):
print(`GWtE(2,2,z,t);`):


elif nops([args])=1 and op(1,[args])=Milim then
print(`Milim(m,k): all the words of length k in the alphabet {1, ...,m}`):
print(`For example, try:`):
print(`Milim(2,4);`):

elif nops([args])=1 and op(1,[args])=Mishpat then
print(`Mishpat(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m,`):
print(`and two words w1 and w2 of the same length in the`):
print(`alphabet {1, ..., m}, and symbols n and N`):
print(`P and t.  It also inputs a positive integer`):
print(`K (the number of terms for guessing the alg. relation, see below).`):
print(`We are interested in the integer sequence`):
print(`a(n):=number of n-lettered words in the alphabet {1, ...,m}`):
print(`that has as many substrings (alias factors, `):
print(`alias consecutive subwords) that are w1 as there are w2.`):
print(`The output this procedure  is a verbose theorem and sketch of proof`):
print(`about the sequence a(n)`):
print(`For example, for Richard Stanley's Amer. Math. Monthly Problem`):
print(`11610 (118(10)[Dec. 2011], p. 937)`):
print(`type: `):
print(`Mishpat(2,[1,1],[1,2],n,N,P,t,40,3);`):

elif nops([args])=1 and op(1,[args])=ProveGF then
print(`ProveGF(m,w1,w2,z,F): Given a positive integer`):
print(`m and two words (lists),w1,w2, of the same length in the`):
print(`alphabet {1, ...,m}, and an expression F that`);
print(`depends on z,`):
print(`proves (rigorously!), using`):
print(`the continuous Almkvist-Zeilberger algorithm,`):
print(`that F is the ordinary generating function (using `):
print(`the variable z) of the sequence`):
print(`"number of words in {1, ...m} of length n that has`):
print(`as many factors (substrings) that are w1 and`):
print(`factors that are w2 equals `):
print(`such that the constant term of s is the generating`):
print(`function for words in {1,..,m}, according to length that`);
print(`have the same number of occurences of the factor w1 and w2.`):
print(`For example, to completely solve Richard Stanley's`):
print(`Amer. Math. Monthly problem #11610 (page 937`):
print(` of v. 118 (No. 10), [Dec. 2011]), please type:`):
print(`ProveGF(2,[1,1],[1,2],t,1/(1-t)/2+(1+2*t)/sqrt((1-t)*(1-2*t)*(1+t+2*t^2))/2);`):

elif nops([args])=1 and op(1,[args])=RECaz then
print(`RECaz(m,w1,w2,n,N): Given a positive integer`):
print(`m and two words (lists),w1,w2, of the same length in the`):
print(`alphabet {1, ...,m}, `):
print(`uses the continuous Almkvist-Zeilberger algorithm,`):
print(`to find a linear recurrence operator with`):
print(`polynomial coefficients, ope(n,N) in terms of n and the`):
print(`shift operator in n, N, that annihilates`):
print(`the sequence`):
print(`"number of words in {1, ...m} of length n that has`):
print(`as many factors (substrings) that are w1 and`):
print(`factors that are w2 equals `):
print(`For example, to get a recurrence operator for`):
print(`Richard Stanley's`):
print(`Amer. Math. Monthly problem #11610 (page 937`):
dprint(` of v. 118 (No. 10), [Dec. 2011]), please type:`):
print(`RECaz(2,[1,1],[1,2],n,N);`):


elif nops([args])=1 and op(1,[args])=RECg then
print(`RECg(m,w1,w2,n,N,MaxC): Given a positive integer`):
print(`m and two words (lists),w1,w2, of the same length in the`):
print(`alphabet {1, ...,m}, `):
print(`uses guessing`):
print(`to find a linear recurrence operator with`):
print(`polynomial coefficients, ope(n,N) in terms of n and the`):
print(`shift operator in n, N`):
print(`such that the degree plus order is <=MaxC ,`):
print(` that annihilates the sequence`):
print(`"number of words in {1, ...m} of length n that has`):
print(`as many factors (substrings) that are w1 and`):
print(`factors that are w2 equals `):
print(`It also gives the initial values.`):
print(`For example, to get a recurrence operator for`):
print(`Richard Stanley's`):
print(`Amer. Math. Monthly problem #11610 (page 937`):
dprint(` of v. 118 (No. 10), [Dec. 2011]), please type:`):
print(`RECg(2,[1,1],[1,2],n,N,10);`):

elif nops([args])=1 and op(1,[args])=RSct then
print(`RSct(m,w1,w2,z,s): The rational function in z and s`):
print(`such that the constant term of s is the generating`):
print(`function for words in {1,..,m}, according to length that`):
print(`have the same number of occurences of the factor w1 and w2.`):
print(`For example, try:`):
print(`RSct(2,[1,1],[1,2],z,s);`):

elif nops([args])=1 and op(1,[args])=RSseq then
print(`RSseq(m,w1,w2,N): Given a pos. integer m and and`):
print(`two words of the same length w1,w2 in the alphabet `):
print(`{1,...,m}, finds the first N terms of the sequence`):
print(`enumerating the number of n-letter words in the`):
print(`alphabet {1,...,m} that have the same number of `):
print(`factors that are w1 and that are factors that are w2`):
print(`For example, try:`):
print(`RSseq(2,[1,1],[1,2],10);`):

 
elif nops([args])=1 and op(1,[args])=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-2,n,N,[1],10);`):

elif nops([args])=1 and op(1,[args])=TerseBook then
print(`TerseBook(m,k,n,N,P,t,K): inputs positive integers m and k`):
print(`and outputs a set of lists consiting of `):
print(`the input and output of Info1(m,k,n,N,P,t,K)`):
print(`(qv) for all members of Zugot(m,k) (qv).`):
print(`For example, try:`):
print(`TerseBook(2,2,n,N,P,t,40);`):

elif nops([args])=1 and op(1,[args])=VerboseBook then
print(`VerboseBook(m,k,n,P,t,K,M): inputs positive integers m and k`):
print(`and outputs a web-book with generating functions`):
print(`and recurrences for sequences enumerating`):
print(`analogous quantities to Richard Stanley's Monthly`):
print(`Problem 11610 (118(10)[Dec. 2011], p. 937)`):
print(`concerning all pairs of words in {1, ...,m} `):
print(`of length k (up to equivalence). `):
print(`m is the size of the alphabet, k is the length of the`):
print(`distinguished words, n,P,t is the symbol for displaying`):
print(`K and M are as in Info2 `):
print(`For example try:`):
print(`VerboseBook(2,2,n,P,t,40,4);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(w,k,z,t): The weight of the word w in terms of z and`):
print(`the indexed variables t by factors of length k. For`):
print(`example, try:`):
print(`Wt([1,2,3,2,1],2,z,t);`):

elif nops([args])=1 and op(1,[args])=WtE then
print(`WtE(m,k,z,t,S): The  weight-enumerator of words in`):
print(`the alphabet {1,..,m} according to the weight z^length`):
print(`times the product of t[FactorsOfLength_kThatAreInTheSetS]. For`):
print(`example, try:`):
print(`WtE(2,2,z,t,{[1,1],[1,2]});`):

elif nops([args])=1 and op(1,[args])=WtEd then
print(`WtEd(m,k,z,t): Like WtE(m,k,z,t,S) `):
print(`but done directly. For checking purposes only. For example, try:`):
print(`WtEd(2,2,z,t,{[1,1],[1,2]});`):

elif nops([args])=1 and op(1,[args])=Zugot then
print(`Zugot(m,k): Representatives of sets of pairs of words of`):
print(`length k in the alphabet {1, ..., m}. For example, try:`):
print(`Zugot(2,2);`):

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





#########Begin AsyRec#########

##From AsyRec
with(SolveTools):
with(numtheory):

#MonoN(ope,n,N,k): Given a linear recurrence operator with
#poly coeffs. ope(n,N), and a pos. integer k, outputs the expression of
#N^k as a linear combination of 1, N, ..., N^(ORDER-1)
#For example, try
#MonoN(N-n-1,n,N,3);
MonoN:=proc(ope,n,N,k) local ORDER,coe0,i,lu1,lu2:
ORDER:=degree(ope,N):

if k<ORDER then
 RETURN(N^k):
fi:

if k=ORDER then
coe0:=coeff(ope,N,ORDER):
RETURN(add(normal(-coeff(ope,N,i)/coe0)*N^i,i=0..ORDER-1)):
fi:

lu1:=expand(MonoN(ope,n,N,k-1)):
lu2:=expand(MonoN(ope,n,N,ORDER)):
lu1:=expand(subs(n=n+1,lu1)*N):
lu1:=subs(N^ORDER=lu2,lu1):
add(normal(coeff(lu1,N,i))*N^i,i=0..ORDER-1):

end:


#Ksect(ope,n,N,k,r): Given a linear recurrence operator ope(n,N)
#annihilating a sequence a[n], and pos. integer k, and
#non-neg. integer r, outputs the one, of the same
#order, that annihilates  a[k*n+r]
#For example try: Ksect(N^2-N-(n+1),n,N,2,0);
Ksect:=proc(ope,n,N,K,r) local ORDER,eq,var,lu,mu,c,Ope,var1,i,vu,v:
 
ORDER:=degree(ope,N):

for i from 0 to ORDER do
 lu[i]:=MonoN(ope,n,N,r+K*i):
 lu[i]:=subs(n=n*K,lu[i]):
od:

Ope:=add(c[i]*N^i,i=0..ORDER):

mu:=expand(add(c[i]*lu[i],i=0..ORDER)):

eq:={seq(coeff(mu,N,i),i=0..ORDER)}:
var:={seq(c[i],i=1..ORDER)}:

var1:=solve(eq,var):


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



vu:={}:

for i from 1 to nops(var1) do
 if op(1,var1[i])=op(2,var1[i]) then
    vu:=vu union {op(1,var1[i])}:
 fi:
od:

Ope:=subs(var1,Ope):

for v in vu do
Ope:=subs(v=0,Ope):
od:



Ope:=subs(c[0]=1,numer(normal(Ope))):



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
#For example, try:
#SeqFromRec(N-2,n,N,[1],10);
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):


if HSP(ope1,n,N)>=1 then
 RETURN(FAIL):
fi:

if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

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

gu:

end:
 



#HafelOper(ope,n,N,L): applies the operator ope(n,N) to
#the sequence L, for example, try:
#HafelOper(N-(n+1),n,N,[1,2,6,24,120]);
HafelOper:=proc(ope,n,N,L) local i,n1:
[seq(add(subs(n=n1,coeff(ope,N,i))*L[n1+i],i=0..degree(ope,N)),
n1=1..nops(L)-degree(ope,N))]:
end:


#TestKsect(ope,n,N,K,r,M,Ini): tests procedure Ksect
#with initial conditions Ini up to M terms
TestKsect:=proc(ope,n,N,K,r,M,Ini)  local gu,Ope,n1:
gu:=SeqFromRec(ope,n,N,Ini,M*K+r):
if gu=FAIL then
  RETURN(FAIL):
fi:
Ope:=Ksect(ope,n,N,K,r):
gu:=[seq(gu[K*n1+r],n1=1..M-1)]:
evalb({op(HafelOper(Ope,n,N,gu))}={0}):
end:

rf:=proc(a,n) local i:mul(a+i,i=0..n-1):end:

#CODV1(ope,n,N,k): Given a linear recurrence operator ope(n,N)
#annihilating a[n], say, outputs the operator
#annihilating b[n]:=a[n]/n!^k
#For example, try: CODV1(N-(n+1),n,N,1):
CODV1:=proc(ope,n,N,k) local i:
add(coeff(ope,N,i)*(rf(n+1,i))^k*N^i,i=0..degree(ope,N)):
end:


#CODV(ope,n,N,k,K): Given a linear recurrence operator ope(n,N)
#annihilating a[n], say, outputs the operator
#annihilating b[n]:=a[n*K]/n!^k
#For example, try: CODV(N-(n+1),n,N,1,1):
CODV:=proc(ope,n,N,k,K) local Ope:
Ope:=Ksect(ope,n,N,K,0):
CODV1(Ope,n,N,k):
end:



#FindKk(ope,n,N): Given a linear recurrence operator with polynomial
#coefficients, ope(n,N), finds the integers K and k such
#that if a(n) is a solution of ope(n,N)a(n)=0, then
#b(n):=a(K*n)/n!^k is annihilated by a standard operator
#the output is the pair [k,K] and the transformed operator
#For exanple, try: FindKk(N^2-n,n,N);
FindKk:=proc(ope,n,N) local ope1,Ld,deg,sp,yakhas,K,k,pu:


pu:=Aope1(ope,n,N):

if {solve(pu,N)}<>{} and {solve(pu,N)}<>{0} then
 RETURN([1,0],ope):
fi:

ope1:=expand(ope):
Ld:=ldegree(ope,N):
ope1:=expand(subs(n=n-Ld,ope1)/N^Ld ):

deg:=degree(ope,N):

sp:=degree(coeff(ope,N,0),n)-degree(coeff(ope,N,deg),n):

yakhas:=sp/deg:

k:=numer(yakhas):
K:=denom(yakhas):
[K,k],CODV(ope,n,N,k,K):
end:


#Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N)
Aope1:=proc(ope,n,N)
local gu:
gu:=expand(numer(normal(ope))):
gu:=coeff(gu,n,degree(gu,n)):
gu:=expand(gu/coeff(gu,N,degree(gu,N))):
factor(gu):
end:

C3:=proc(a,b,c) option remember:

if [a,b,c]=[0,0,0] then
  RETURN(1):
fi:
if a<0 or b<0 or c<0 then
   RETURN(0):
fi:
C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1):
end:

#Aluf1(pit): is it positive dominant?
Aluf1:=proc(pit)
local aluf,shi,i,makom:


shi:=evalf(abs(pit[1])):
aluf:={evalf(evalc(pit[1]))}:
makom:={1}:

for i from 2 to nops(pit) do

 if evalf(abs(pit[i]))>evalf(shi) then
   shi:=abs(evalf(evalc(pit[i]))):
  aluf:={pit[i]}:
  makom:={i}:
 elif evalf(abs(pit[i]))=shi then
   aluf:=aluf union {evalf(evalc(pit[i]))}:
   makom:=makom union {i}:
 fi:

od:

aluf,shi,makom:

end:



#OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic
#expansion of solutions of ope1(n,N) with exponent alpha
#extends it by one term
OneStepAS1:=proc(ope1,n,N,alpha,f,S1)
local x,f1,L,F,A,mu,i,A1:

f1:=subs(n=1/x,f):
L:=degree(f1,x):
F:=f1+A*x^(L+S1):
mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F)
,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
mu:=taylor(mu,x=0,L+11):
mu:=simplify(mu):
for i from L+1 to L+9 while coeff(mu,x,i)=0 do od:
if i=L+10 then
 RETURN(FAIL):
fi:
mu:=coeff(mu,x,i):
A1:=subs(Linear({mu},{A}),A):
if A1=0 then
  RETURN(FAIL):
fi:


subs({x=1/n,A=A1},F):

end:


#OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic
#expansion of solutions of ope1(n,N) with exponent alpha
#extends it by one term
OneStepA:=proc(ope1,n,N,alpha,f)
local S1,pu:
for S1 from 1 to 5 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od:
pu:=OneStepAS1(ope1,n,N,alpha,f,S1):
if pu=FAIL then
 RETURN(f):
else
 RETURN(pu):
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:






#Nor(ope,n,N): the Normalizer of the linear recurrence
#operator with polynomial coefficients ope(n,N)
#followed by its exponential growth constant
#For example, try:
#Nor(N^2-3*N+2,n,N);
Nor:=proc(ope,n,N)
local gu,pit,lu,makom,ope1:

gu:=Aope1(ope,n,N):
if degree(gu,N)=0 then
   print(`Not of exponential growth, case not implemented`):
   RETURN(FAIL):
fi:

if not type(expand(normal(ope/gu)),`+`) then
 pit:=[solve(gu,N)]:
 makom:=Aluf1(pit)[3]:

lu:=pit[makom[1]]:
ope1:=Yafe(subs(N=lu*N,ope),N)[2]:
RETURN(ope1,lu):

fi:

pit:=[solve(gu,N)]:
makom:=Aluf1(pit)[3]:

if nops(makom)<>1 then
lu:=evalf(evalc(pit[makom[1]])):
 if coeff(lu,I,1)>(1/10)^(Digits-2) then
 RETURN(FAIL):
  elif pit[makom[1]]+pit[makom[2]]=0 then
   RETURN(FAIL):
 fi:
fi:


lu:=evalf(evalc(pit[makom[1]])):

if coeff(lu,I,1)>(1/10)^(Digits-2) then
 RETURN(FAIL):
fi:

lu:=pit[makom[1]]:

ope1:=Yafe(subs(N=lu*N,ope),N)[2]:
ope1,lu:
end:









#Atom(s,r,i,x,K): Expanding (n+i)^(s/r)-n^(s/r) in terms of
#x=n^(-1/r) up to K terms
#For example try:
#Atom(2,3,2,x,3);
Atom:=proc(s,r,i,x,K) local gu,y,i1:
gu:=(1+i*y)^(s/r)-1:
gu:=taylor(gu,y=0,K+1):
gu:=add(coeff(gu,y,i1)*y^i1,i1=0..K):
expand(subs(y=x^r,gu)/x^s):
end:




#FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics
#in terms of x=n^(1/r)
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x);
FindExpP1old:=proc(ope,n,N,r,x) local c,eq,var,s,sof,gu,ope1,deg,i,i1,v:
sof:=add(c[s]*x^s,s=1..r-1):
var:={seq(c[i],i=1..r-1)}:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=expand(subs(n=1/x^r,ope1)):

gu:=0:

for i from 0 to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..r-1) ):
od:

gu:=taylor(gu,x=0,2*r-1):



eq:={seq(coeff(gu,x,i1),i1=r..2*r-2)}:


var:=[solve(eq,var)][1]:
for v  in var do
 if op(1,v)=op(2,v) then
  RETURN(FAIL):
 fi:
od:

subs(var,sof):

end:







#Findr(ope,n,N): Given an operator ope(n,N)
#such that the leading terms in n is a multiple of
#(N-1), finds the largest r such that it is
# a multiple of (N-1)^r. For example, try:
#Findr((N-1)^4,n,N);
Findr:=proc(ope,n,N) local ope1,r:
ope1:=coeff(expand(ope),n,degree(ope,n)):
if expand(subs(N=1,ope1))<>0 then
  RETURN(FAIL):
fi:

for r from 1 while expand(subs(N=1,diff(ope1,N$r)))=0 do od:
r:
end:


FindExpP:=proc(ope,n,N)  local r,x,lu:
 r:=Findr(ope,n,N):


  if r=FAIL then
    RETURN(FAIL):
  fi:
 lu:=FindExpP1(ope,n,N,r,x) :
  if lu=FAIL then
    RETURN(FAIL):
   fi:
  subs(x=n^(1/r),lu):
end:




#NewOpe(ope,n,N,K): the exponential part + the transformed operator
#(up to degree K asymp. in the coefficients)
#For example, try:
#NewOpe((n+1)*N^2-2*(n+5)*N+n+3,n,N,4);
NewOpe:=proc(ope,n,N,K,x)  local r,lu,lu1,ope1,deg,s,i,gu,mu:
 r:=Findr(ope,n,N):
  if r=FAIL then
    RETURN(FAIL):
  fi:
 lu:=FindExpP1(ope,n,N,r,x) :
 lu1:=FindExpP(ope,n,N) :

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

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=subs(n=1/x^r,ope1):
ope1:=expand(ope1):


gu:=0:

for i from 0 to degree(ope1,N) do
  mu:=0:
 for s from 1 to r-1 do
     mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3):
  od:
     
 gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i:
od:

gu:=taylor(gu,x=0,K+3):
lu1,add(coeff(gu,x,i)*x^i,i=0..K):
end:



#Bdok(ope,n,N,K): the exponential part + the transformed operator
#(up to degree K asymp. in the coefficients)
#For example, try:
#Bdok((n+1)*N^2-2*(n+5)*N+n+3,n,N,4);
Bdok:=proc(ope,n,N,K)  local r,x,lu,lu1,ope1,deg,s,i,gu,mu:
 r:=Findr(ope,n,N):
  if r=FAIL then
    RETURN(FAIL):
  fi:
 lu:=FindExpP1(ope,n,N,r,x) :
 lu1:=FindExpP(ope,n,N) :
  if lu=FAIL then
    RETURN(FAIL):
   fi:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=subs(n=1/x^r,ope1):
ope1:=expand(ope1):


gu:=0:

for i from 0 to degree(ope1,N) do
  mu:=0:
 for s from 1 to r-1 do
     mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3):
  od:
     
 gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i:
od:

gu:=taylor(gu,x=0,K+3):
lu1,add(coeff(gu,x,i)*x^i,i=0..K):
end:




OpePerN:=proc(N,n,r) local i,j,ope:
ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r):
ope:=expand(N^r*subs(n=n+r,ope)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r):
ope:=expand(FindKk(ope,n,N)[2]):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r):
ope:=numer(Nor(ope,n,N)[1]):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r):
end:

OpePer:=proc(N,n,r) local i,j,ope:
ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r):
ope:=expand(N^r*subs(n=n+r,ope)):
end:

OpePerG:=proc(N,n,S) local i,j,ope,r,gu,t:
ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i in S):
r:=max(op(S)):
ope:=expand(N^r*subs(n=n+r,ope)):

gu:=exp(add(t^i/i,i in S)):
gu:=taylor(gu,t=0,degree(ope,N)+2):

ope,[seq(coeff(gu,t,i)*i!,i=1..degree(ope,N))]:
end:


#Finda(ope,N,x,r): finds the first power x^a in the
#asymptotic solution of ope(N,n)f(n)=0, where x=1/n^(1/r)
#For example, try:
#Finda((1+x)-(1+3*x)*N,N,x,1);
Finda:=proc(ope,N,x,r) local gu,i,a,a1:

gu:=add(coeff(ope,N,i)*(1+i*x^r)^(a/r),i=0..degree(ope,N)):

gu:=taylor(gu,x=0,5*r+5):
gu:=expand(add(coeff(gu,x,i)*x^i,i=0..5*r+4)):

for i from 0 to 5*r+4 while expand(simplify(coeff(gu,x,i)))=0 do od:

if i=5*r+5 then
  RETURN(FAIL):
fi:


a1:=[solve(coeff(gu,x,i),a)][1]:


if a1=NULL then
 RETURN(FAIL):
elif not type(a1,integer) then
 RETURN(FAIL):

else
 RETURN(-a1):
fi:
end:






#OneStepG(ope,N,x,r,Cu): Given a partial
#asympt. expansion, in terms of x=1/n^(1/r),
#for solutions of the recurrence equation
#ope(n,N)a(n)=0, where ope is normalized of type
#r (i.e. its leading coeff. is a multiple of (N-1)^r)
#finds the next term
#OneStepG((1+x)-(1+3*x)*N,N,x,1,x^2);
OneStepG:=proc(ope,N,x,r,Cu) local gu,i,a,c,c1,Cu1,K,ka:
K:=degree(Cu,x):
Cu1:=Cu+c*x^(K+1):
gu:=add(coeff(ope,N,i)*
add(coeff(Cu1,x,a)*x^a*(1+i*x^r)^(-a/r),a=ldegree(Cu1,x)..degree(Cu1,x)),
i=0..degree(ope,N)):
gu:=expand(gu):

gu:=taylor(gu,x=0,5*r+8+K):
gu:=
simplify(expand(add(coeff(gu,x,i)*x^i,i=0..5*r+7+K))):

gu:=pashet(gu,x,5*r+7+K,c):
ka:=expand(simplify(coeff(gu,x,ldegree(gu,x)))):
ka:=simplify(ka):

c1:=[solve(simplify(coeff(gu,x,ldegree(gu,x))),c)][1]:


if c1=NULL then
 RETURN(FAIL):
else
 RETURN(subs(c=c1,Cu1)):
fi:
end:



#Asy1(ope,n,N,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the K's term
Asy1:=proc(ope,n,N,K)
local gu,lu,alpha,mu,ope1,ku,i,f,x,ka,vu,ope1A,r,Kk,pu,a:


gu:=FindKk(ope,n,N):
Kk:=gu[1]:

ope1:=gu[2]:

vu:=Nor(ope1,n,N):

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

ope1:=vu[1]:
lu:=vu[2]:

ope1A:=Aope1(ope1,n,N):



for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od:


if r=1 then

###r=1 case

mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
ku:=factor(coeff(taylor(mu,x=0,2),x,1)):


ka:=simplify([solve(ku,alpha)]):


if coeff(ka[1],I,1)<>0 then
  RETURN(FAIL):
fi:

alpha:=max(op(ka)):


if normal(subs(N=1,ope1))=0 then
 RETURN(lu^n*n^alpha,Kk):
fi:

f:=1:
for i from 1 to K do
f:=OneStepA(ope1,n,N,alpha,f):
od:

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

RETURN(lu^n*n^alpha*f,Kk):

#end r=1 case

else



ope1:=numer(ope1):


if expand(diff(ope1,n))=0 then
 ka:=n^(r-1)*lu^n:
 RETURN(ka,Kk):
fi:

 pu:=FindExpP(ope1,n,N):

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

ope1:=NewOpe(ope1,n,N,K+5,x)[2]:

a:=Finda(ope1,N,x,r):

if a=FAIL then
 RETURN(FAIL):
fi:
ka:=x^a:

for i from a to K do
 ka:=OneStepG(ope1,N,x,r,ka):
od:

ka:=exp(pu)*subs(x=1/n^(1/r),ka)*lu^n:

RETURN(ka,Kk):
fi:
end:



pashet:=proc(gu,x,K,c) local gu1,i,pu:

Digits:=300:

gu1:=0:

for i from 0 to K do
 pu:=evalf(coeff(gu,x,i)):
if not (abs(coeff(pu,c,1))<10^(-20)  and abs(coeff(pu,c,0))<10^(-20) )
then
  gu1:=gu1+coeff(gu,x,i)*x^i:
fi:
od:
gu1:
end:


  
 
#NakedStirling(n,K): the asymptotic expansion of
#n!/((n/e)^n*sqrt(2*Pi*n). For example, try:
#NakedStirling(n,5);
NakedStirling:=proc(n,K) local gu,i:
gu:=n!/(n/exp(1))^n/sqrt(2*Pi*n):
gu:=expand(asympt(gu,n,K+1)):

(n/exp(1))^n*sqrt(2*Pi*n), add(coeff(op(i,gu),n,-(i-1))/n^(i-1),i=1..K+1):
end:


#AsyF(ope,n,N,M): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the M's term, in terms of the factorial function
AsyF:=proc(ope,n,N,M) local K,k,gu:
gu:=Asy1(ope,n,N,M):

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

K:=gu[2][1]:
k:=gu[2][2]:
gu:=gu[1]:
(n/K)!^k*subs(n=n/K,gu):
end:



#AsyFC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0, with the given initial conditions
#where ope(n,N) is  a recurrence operator
#up to the M's term, in terms of the factorial function
#and complete with an empirically derived constant in front
#using K terms
AsyFC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1:
Digits:=100:
gu:=AsyF(ope,n,N,M):
L:=SeqFromRec(ope,n,N,evalf(Ini),K):
if L=FAIL then
 RETURN(FAIL):
fi:

mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]:
er:=mu[nops(mu)]-mu[nops(mu)-1]:
D1:=-trunc(log(abs(er))/log(10))-3:
if D1<2 then
 print(`can't determine the constant`):
RETURN(gu):
fi:

C:=evalf(mu[nops(mu)],D1):
C,gu:
end:

#Asy1special(ope,n,N,K,x): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,given as a list
#[pu,lu,expansion,r] where it is
#exp(pu)*lu^n*expansion(x), where x=1/n^(1/r), and r is
#a positive integer. It also returns [K,k] (see Asy1)
#where ope(n,N) is  a recurrence operator
#up to the K's term
Asy1special:=proc(ope,n,N,K,x)
local gu,lu,alpha,mu,ope1,ku,i,f,ka,vu,ope1A,r,Kk,pu,a:


gu:=FindKk(ope,n,N):



Kk:=gu[1]:

ope1:=gu[2]:

vu:=Nor(ope1,n,N):


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

ope1:=vu[1]:
lu:=vu[2]:

ope1A:=Aope1(ope1,n,N):



for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od:


if degree(ope1,n)=0 then
RETURN([0,lu,x^(-r+1),1,1],Kk):
fi:


if r=1 then

###r=1 case

mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
ku:=factor(coeff(taylor(mu,x=0,2),x,1)):


ka:=simplify([solve(ku,alpha)]):

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

if coeff(ka[1],I,1)<>0 then
  RETURN(FAIL):
fi:

alpha:=max(op(ka)):


if normal(subs(N=1,ope1))=0 then
RETURN([0,lu,x^(-alpha),1,1],Kk):
fi:

f:=1:
for i from 1 to K do
f:=OneStepA(ope1,n,N,alpha,f):
od:

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

f:=subs(n=1/x,f):
RETURN([0,lu,x^(-alpha),f,1],Kk):
#end r=1 case

else



ope1:=numer(ope1):


if expand(diff(ope1,n))=0 then
 ka:=n^(r-1)*lu^n:
 RETURN(ka,Kk):
fi:




 pu:=FindExpP(ope1,n,N):



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


ope1:=NewOpe(ope1,n,N,K+5,x)[2]:

a:=Finda(ope1,N,x,r):

if a=FAIL then
 RETURN(FAIL):
fi:
ka:=x^a:

for i from a to K do
 ka:=OneStepG(ope1,N,x,r,ka):
od:


RETURN([pu,lu,1,ka,r],Kk):

fi:
end:

#IsMispar(a): is a (generalized) numeric type
IsMispar:=proc(a) :
if type(a,numeric) then
 true:
elif type(a,`^`) and type(op(1,a),numeric) and type(op(2,a),numeric)  then
 true:
elif a=Pi  then
true:
elif type(a,`^`) and op(1,a)=Pi and type(op(2,a),numeric)  then
 true:
elif type(a,function) and type(op(1,a),numeric) then
 true:
else
 false:
fi:
end:

#GetRidOfConst(P): gets rid of the constant in front
GetRidOfConst:=proc(P) local i,P1:
if not type(P,`*`) then
 RETURN(P):
fi:

P1:=1:

for i from 1 to nops(P) do
 if not IsMispar(op(i,P)) then
  P1:=P1*op(i,P):
 fi:
od:
P1:
end:

#Asy(ope,n,N,M): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the M's term
Asy:=proc(ope,n,N,M) local K,k,gu,vu,vu1,vu2,pu,lu,ka1,ka2,r,x,i,vu2a:
gu:=Asy1special(ope,n,N,M,x):
if gu=FAIL then
 RETURN(FAIL):
fi:

K:=gu[2][1]:
k:=gu[2][2]:


gu:=gu[1]:
pu:=subs(n=n/K,gu[1]):
lu:=gu[2]:
ka1:=gu[3]:
ka2:=gu[4]:
r:=gu[5]:
ka2:=subs(x=K^(1/r)*x,ka2):
vu:=NakedStirling(n,M+1):
vu1:=vu[1]:
vu2:=vu[2]:


vu1:=subs(n=n/K,vu1)^k:
vu2:=subs(n=n/K,vu2)^k:
vu2:=expand(subs(n=1/x^r,vu2)):

vu2:=expand(ka2*vu2):


vu2:=add(simplify(coeff(vu2,x,i))*x^i,i=0..M*r):
vu2:=subs(x=1/n^(1/r),vu2):


ka1:=subs(x=1/n^(1/r),ka1):
ka1:=subs(n=n/K,ka1):
vu2a:=op(1,vu2):
vu2:=expand(normal(vu2/vu2a)):
gu:=simplify(vu2a*vu1*lu^(n/K))*exp(pu)*ka1*vu2:

GetRidOfConst(gu):

end:






#AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0, with the given initial conditions
#where ope(n,N) is  a recurrence operator
#up to the M's term,
#and complete with an empirically derived constant in front
#using K terms
AsyC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1,gu1:
Digits:=100:
gu:=Asy(ope,n,N,M):
if gu=FAIL then
 RETURN(FAIL):
fi:
L:=SeqFromRec(ope,n,N,evalf(Ini),K):
if L=FAIL then
 RETURN(FAIL):
fi:

mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]:
er:=mu[nops(mu)]-mu[nops(mu)-1]:


if er=0 then
 D1:=Digits:
else
D1:=-trunc(evalf(log(abs(er))/log(10)))-3:
fi:

if D1<2 then
 print(`can't determine the constant`):
RETURN(FAIL):
fi:

C:=evalf(mu[nops(mu)],D1):

if C<1/10000 then
 RETURN(FAIL):
fi:

gu1:=C*gu:


if not CheckAsy(gu1,n,L,1/100) then
 RETURN(FAIL):
fi:

C,gu,D1:

end:








#FindExpP1g(ope,n,N,r,x,ds): finds the exponential part of the asymptotics
#in terms of x=n^(1/r) as a poly. of degree ds in x
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1g((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1);
FindExpP1g:=proc(ope,n,N,r,x,ds) local c,eq,var,s,sof,gu,ope1,
deg,i,i1,v,ka,i11:
sof:=add(c[s]*x^s,s=1..ds):
var:={seq(c[i],i=1..ds)}:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=expand(subs(n=1/x^r,ope1)):


gu:=0:

for i from 0 to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ):
od:


gu:=taylor(gu,x=0,3*r+10):


for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od:

eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1)}:
eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds)}:

lprint(eq,var):
ka:=[solve(eq,var)]:


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


var:=ka[1]:
for v  in var do
 if op(1,v)=op(2,v) then
  RETURN(FAIL):
 fi:
od:

subs(var,sof):

end:


#FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics
#in terms of x=n^(1/r)
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x);
FindExpP1:=proc(ope,n,N,r,x) local s,gu,nosaf:

for s from 1 to r-1 do
 for nosaf from 0 to 1 do
 gu:=FindExpP1gExtra(ope,n,N,r,x,s,nosaf):
 if gu<>FAIL then
  RETURN(gu):
 fi:
od:
od:

FAIL:
end:





#FindExpP1gExtra(ope,n,N,r,x,ds,nosaf)
#: finds the exponential part of the asymptotics
#in terms of x=n^(1/r) as a poly. of degree ds in x
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1gExtra((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1);
FindExpP1gExtra:=proc(ope,n,N,r,x,ds,nosaf) local c,eq,var,s,sof,gu,ope1,
deg,i,i1,v,ka,i11,ku,varf:
sof:=add(c[s]*x^s,s=1..ds):
var:={seq(c[i],i=1..ds)}:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=expand(subs(n=1/x^r,ope1)):


gu:=0:

for i from 0 to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ):
od:


gu:=taylor(gu,x=0,3*r+10):


for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od:

eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1+nosaf)}:


ka:=[solve(eq,var)]:


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


var:=ka[1]:
for v  in var do
 if op(1,v)=op(2,v) then
  RETURN(FAIL):
 fi:
od:

#varf:=evalf(var):

#ku:=[seq(subs(varf,c[s]),s=1..ds)]:

#print(`ku is`, ku):
#add(ku[s]*x^s,s=1..ds):

subs(var,sof):

end:
#########End from AsyRec#########
 



###Begin from SCHUTZENBERGER#####

simp:=proc(pa,deg,bitui,P,x)
local mone,mekh,gu:
gu:=normal(bitui):
mone:=numer(gu):
mekh:=denom(gu):
mone:=simp1(pa,deg,mone,P,x):
mekh:=simp1(pa,deg,mekh,P,x):
expand(mone)/expand(mekh):
end:
 
simp1:=proc(pa,deg,bitui,P,x)
local i,gu:
 
gu:=expand(bitui):
 
for i from deg to 3*deg+1 do
 gu:=subs(P(x)^i=pa[i],gu):
od:
 
gu:
 
end:
 

#algtodiff(F,P,x,D1): Inputs a polynomial F,
#in P and x, and outputs
#a differential operator annihilating
#P(x) where P(x) is the alg. formal power series
#satisfying the algebraic equation
#F(P(x),x))=0 . For example, try:
#algtodiff(P-1-x*P^2,P,x);
algtodiff:=proc(F,P,x,D1) local ORDER,gu:

for ORDER from 1 do

gu:=algtodiff1(F,P,x,D1,ORDER):
if gu<>FAIL then
gu:=add(factor(coeff(gu,D1,i))*D1^i,i=0..ORDER):
 RETURN(gu):
fi:
od:

end:

#algtodiff1(F,P,x,D1,ORDER): Inputs a polynomial F,
#in P and x, and an ORDER, and outputs
#a differential operator annihilating
#P(x) where P(x) is the alg. formal power series
#satisfying the algebraic equation
#F(P(x),x))=0 . For example, try:
#algtodiff1(P-1-x*P^2,P,x,2);
algtodiff1:=proc(F,P,x,D1,ORDER)
local i,gu,a,lu,pa,deg,eq1,lu1,y,var,eq,ope,pip:
deg:=degree(F,P):
pa:=array(deg..3*deg+1):
eq1:=subs(P^deg=y,F):
pa[deg]:=solve(eq1=0,y):
pa[deg]:=subs(P=P(x),pa[deg]):
 
for i from 1 to 2*deg+1 do
   lu:=pa[deg+i-1]*P(x):
   lu:=expand(lu):
   lu:=subs(P(x)^deg=pa[deg],lu):
   pa[deg+i]:=lu:
od:
 
 
gu:=a[0]*P(x):
lu:=subs(P=P(x),F):
lu1:=diff(lu,x):
lu1:=subs(diff(P(x),x)=y,lu1):
lu1:=solve(lu1=0,y):
lu1:=simp(pa,deg,lu1,P,x):
 lu1:=normal(lu1):
 
 
gu:=gu+a[1]*lu1:
gu:=normal(gu):
 gu:=simp(pa,deg,gu,P,x):
 
var:={a[0],a[1]}:
ope:=a[0]+a[1]*D1:
lu:=lu1:
 for i from 2 to ORDER do
 lu:=diff(lu,x):
  lu:=subs(diff(P(x),x)=lu1,lu):
  lu:=normal(lu):
 lu:=simp(pa,deg,lu,P,x):
 gu:=gu+a[i]*lu:
 gu:=normal(gu):
 gu:=simp(pa,deg,gu,P,x):
 var:=var union {a[i]}:
 
 ope:=ope+a[i]*D1^i:
od:
 
gu:=normal(gu):
 
 
gu:=numer(gu):
 
gu:=expand(gu):
eq:={}:
 
for i from 0 to deg-1 do
 eq:=eq union {coeff( gu,P(x),i)=0}:
od:
 
 
var:=solve(eq,var):
ope:=subs(var,ope):
ope:=normal(ope):
ope:=numer(ope):

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

ope:=expand(ope):
pip:=expand(coeff(ope,D1,degree(ope,D1))):

pip:=coeff(pip,x,degree(pip,x)):


normal(ope/pip):
end:

#empir(gu,degx,degP,x,P) 
#to "fit" an algebraic equation F(P(x),x)=0 of degree
#degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i
 
empir:=proc(gu,degx,degP,x,P)
local i1,i2,F,a,cand,lu,eq,var,mu,ka,mu1,K,richard,ka1:
K:=(1+degx)*(1+degP)+6:
if K > nops(gu) then
 print(`sequence too small`):
 RETURN(FAIL):
fi:
 
F:=0:
var:={}:
for i1 from 0 to degx do
 for i2 from 0 to degP do
  F:=F+a[i1,i2]*x^i1*P^i2:
  var:=var union {a[i1,i2]}:
 od:
od:
 
cand:=0:
 
for i1 from 0 to nops(gu)-1 do
 cand:=cand+op(i1+1,gu)*x^i1:
od:
 
lu:=subs(P=cand,F):
lu:=taylor(lu,x=0,K+1):
 
eq:={}:
 
for i1 from 0 to K do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):

ka:={}:

for mu1 in mu do
 if op(1,mu1)=op(2,mu1) then
  ka:=ka union {op(1,mu1)}:
 fi:
od:

 
F:=subs(mu,F):
 
if F=0 then
 RETURN(0):
fi:

richard:={seq(factor(coeff(F,ka1,1)) ,ka1 in ka)}:

if nops(richard)>1 then
RETURN(0):
else
 RETURN(richard[1]):
fi:

end:

#Empir(L,x,P): inputs a sequence L and 
#and conjectures an polynomial equation
#F(x,P)=0 satisfies by the formal power
#series L[1]+L[2]*x+L[3]*x^2+...
#For example, try:
#Empir([seq(i1,i1=1..40)],x,P);
Empir:=proc(L,x,P)
local degx,degP,lu,L1,i:
 
for L1 from 1 to (nops(L)-3)/3 do
for degP from 1 to L1 do
for degx from 0 to min(trunc((nops(L)-6)/(1+degP))-1,L1-degP) do
 
lu:=empir(L,degx,degP,x,P):
 
if lu<>0 then
RETURN(add(factor(coeff(lu,P,i))*P^i,i=0..degP)):
fi:
od:
od:
od:
0:
end:

###End from SCHUTZENBERGER#####

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

ope:=subs(var1,ope):

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

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

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

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a,ope1:
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(FAIL):
elif nops(ope)=1 then

if HSP(ope[1],n,N)>=1 then
 RETURN(FAIL):
else
ope1:=Yafe(ope[1],N)[2]:

if not CheckRec(f,ope1,n,N) then
 RETURN(FAIL):
else
RETURN(ope1):
fi:

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

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

end:





 
#SeqFromRecOld(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
SeqFromRecOld:=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)]:
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#########

###Begin From EKHAD
pashetc:=proc(p,D)
local gu,p1,i,mekh:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,D) do
gu:=gu+factor(coeff(p1,D,i))*D^i:
od:
RETURN(gu,mekh):
end:

gzor:=proc(f,x,n)
local i,gu:
gu:=f:
for i from 1 to n do
gu:=diff(gu,x):
od:
RETURN(gu):
end:


goremc:=proc(D,ORDER)
local i,gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*D^i:
od:
RETURN(gu):
end:

duisc:= proc(A,ORDER,y,x,D)
local KAMA,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,i,shad:
 
KAMA:=50:
gorem:=goremc(D,ORDER):
 
conj:=gorem:
yakhas:=0:
 
for i from 0 to degree(conj,D) do
yakhas:=yakhas+normal(coeff(conj,D,i)*gzor(A,x,i)/A):
yakhas:=normal(yakhas):
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:=pashetc(gorem,D):
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:
shad[1],S:
end:


Mygcd:=proc(L) local i,L1:
if nops(L)=1 then
 RETURN(1):
fi:
if nops(L)=2 then
 RETURN(gcd(L[1],L[2])):
fi:
L1:=[op(1..nops(L)-1,L)]:
gcd(Mygcd(L1),L[nops(L)]):
end:

AZc:=proc(A,y,x,D)
local ORDER,gu,KAMA,ope,C,kak,i1:
KAMA:=50:
 
for ORDER from 1 to KAMA do
gu:=duisc(A,ORDER,y,x,D):


if gu[1]<>0 then
ope:=gu[1]:
C:=gu[2]:
kak:=Mygcd([seq(coeff(ope,D,i1),i1=0..ORDER)]):
if kak<>1 then
ope:=add(normal(coeff(ope,D,i1)/kak)*D^i1,i1=0..ORDER):
C:=normal(C/kak):
fi:

 RETURN(ope,C):
fi:
od:
0:
end:

###End From EKHAD

####Start DiffToREc1
#DiffToRec1(ope,t,Dt,n,N): Given a linear differential operator
#ope(t,Dt) and symbols n,N, outputs the linear recurrence
#operator with polynomial coefficients satisfied by
#the sequence of coeff. of any f.p.s. annihilating ope(t,Dt)
#where N is the shift in n. For example, try:
#DiffToRec1(Dt+t, t,Dt,n,N);
DiffToRec1:=proc(ope,t,Dt,n,N) local gu,i,j,ope1,d,j1:
ope1:=expand(ope):

gu:=0:

for i from ldegree(ope1,t) to degree(ope1,t) do
 for j from 0 to degree(ope1,Dt) do
   gu:=gu+coeff(coeff(ope1,t,i),Dt,j)*mul(n-i+j1,j1=1..j)*N^(j-i):
 od:
od:


d:=ldegree(gu,N):
gu:=N^(-d)*subs(n=n-d,gu):



Yafe(gu,N)[2]:

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:

####End DiffToRec1

#Milim(m,k): all the words of length k in the alphabet {1, ...,m}
#For example, try:
#Milim(2,4);
Milim:=proc(m,k) local mu,i,w:
if k=0 then
 RETURN({[]}):
fi:
mu:=Milim(m,k-1):

{seq(seq([op(w),i],i=1..m),w in mu)}:

end:


#EqsVars(m,k,z,t,A):The set of equations, followed by the
#set of variables that enables to compute the generating
#function in z and t[w1] (w1 in [1,...,m]^k)
#that is the weight-enumerator of all words in
#the alpbhabet {1, ...,m} with weight z^length(w) times
#product of t[w_i,...,w_i+k-1] for i from 1 to nops(w)-k+1
#A is a symbol for the indexing of the weight enumerator
#of words,
#For example, try
#EqsVars(2,2,z,t,A);
EqsVars:=proc(m,k,z,t,A) local eq,var,gu,eq1,i,w:
gu:=Milim(m,k-1):
var:={}:
eq:={}:

for w in gu do
 var:=var union {A[op(w)]}:
 eq1:=A[op(w)]-z^(k-1)-z*add(t[op(w),i]*A[op(2..nops(w),w),i],i=1..m):
 eq:=eq union {eq1}:
od:

eq,var:

end:


#GWtE(m,k,z,t): The global weight-enumerator of words in
#the alphabet {1,..,m} according to the weight z^length
#times the product of t[FactorsOfLength_k]. For
#example, try:
#GWtE(2,2,z,t);
GWtE:=proc(m,k,z,t) local gu,A,mu,lu,i1,mu1:
 mu:=EqsVars(m,k,z,t,A):
 gu:=solve(mu):
lu:=add(subs(gu,mu1),mu1 in mu[2]):
lu:=lu+add(m^i1*z^i1,i1=0..k-2):
normal(lu):
end:


#Wt(w,k,z,t): The weight of the word w in terms of z and
#the indexed variables t by factors of length k. For
#example, try:
#Wt([1,2,3,2,1],2,z,t);
Wt:=proc(w,k,z,t) local i:
z^nops(w)*mul(t[op(i..i+k-1,w)],i=1..nops(w)-k+1):
end:




#CheckGWtE(m,k,N): checks the first N terms of
#the output of CheckGWtE(m,k,z,t); For example, try:
#CheckGWtE(2,2,5);
CheckGWtE:=proc(m,k,N) local z,t,gu,mu,lu,w,i:
gu:=GWtE(m,k,z,t):
gu:=taylor(gu,z=0,N+1):

for i from 0 to N do
lu:=expand(coeff(gu,z,i))*z^i:
mu:=Milim(m,i):
if expand(add(Wt(w,k,z,t),w in mu)-lu)<>0 then
 RETURN(false):
fi:
od:
true:

end:


#WtE(m,k,z,t,S): The  weight-enumerator of words in
#the alphabet {1,..,m} according to the weight z^length
#times the product of t[FactorsOfLength_kThatAreInTheSetS]. For
#example, try:
#WtE(2,2,z,t,{[1,1],[1,2]});
WtE:=proc(m,k,z,t,S) local gu,A,mu,lu,i1,mu1,S1,eq,var,s1:
 mu:=EqsVars(m,k,z,t,A):
eq:=mu[1]:
var:=mu[2]:
S1:=Milim(m,k) minus S:
eq:=subs({seq(t[op(s1)]=1,s1 in S1)},eq):
 gu:=solve(eq,var):
lu:=add(subs(gu,mu1),mu1 in mu[2]):
lu:=lu+add(m^i1*z^i1,i1=0..k-2):
normal(lu):
end:


#WtEd(m,k,z,t,S): Like WtE(2,2,z,t,{[1,1],[1,2]});
#but doing it directly, for checking purposes only
#Do
#WtEd(2,2,z,t,{[1,1],[1,2]});
WtEd:=proc(m,k,z,t,S) local gu,S1,s1:
S1:=Milim(m,k) minus S:
gu:=GWtE(m,k,z,t):
gu:=subs({seq(t[op(s1)]=1,s1 in S1)},gu):
normal(gu):
end:

#RSseq(m,w1,w2,N): Given a pos. integer m and and
#two words of the same length w1,w2 in the alphabet 
#{1,...,m}, finds the first N terms of the sequence
#enumerating the number of n-letter words in the
#alphabet {1,...,m} that have the same number of 
#factors that are w1 and that are factors that are w2
#For example, try:
#RSseq(2,[1,1],[1,2],10);
RSseq:=proc(m,w1,w2,N) local z,t,gu,i1,s,L,i,ku:

if not type(m, integer) and m>=1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if not type(N, integer) and N>=1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if not (type(w1,list) and type(w2,list)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(w1)<>nops(w2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {op(w1)} minus {seq(i1,i1=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {op(w2)} minus {seq(i1,i1=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:


gu:=WtE(m,nops(w1),z,t,{w1,w2}):

gu:=subs({t[op(w1)]=s,t[op(w2)]=1/s},gu):
gu:=taylor(gu,z=0,N+1):

L:=[]:
for i from 1 to N do
ku:=expand(coeff(gu,z,i)):
L:=[op(L),coeff(ku,s,0)]:
od:
L:

end:




#RSct(m,w1,w2,z,s): The rational function in z and s
#such that the constant term of s is the generating
#function for words in {1,..,m}, according to length that
#have the same number of occurences of the factor w1
#and w2.
#For example, try:
#RSct(2,[1,1],[1,2],z,s);
RSct:=proc(m,w1,w2,z,s) local gu,t,i1:


if not type(m, integer) and m>=1 then
 print(`Bad input`):
 RETURN(FAIL):
fi:


if not (type(w1,list) and type(w2,list)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(w1)<>nops(w2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {op(w1)} minus {seq(i1,i1=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if {op(w2)} minus {seq(i1,i1=1..m)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:


if not (type(w1,list) and type(w2,list) and 
{op(w1)} minus {seq(i1,i1=1..m)}={} and
{op(w2)} minus {seq(i1,i1=1..m)}={} and
w1<>w2) then
 print(`bad input`):
 RETURN(FAIL):
fi:

gu:=WtE(m,nops(w1),z,t,{w1,w2}):

normal(subs({t[op(w1)]=s,t[op(w2)]=1/s},gu)):

end:




#DOaz(m,w1,w2,z,Dz): Given a positive integer
#m and two words (lists),w1,w2, of the same length in the
#alphabet {1, ...,m}, 
#uses the continuous Almkvist-Zeilberger algorithm,
#to find a linear differential operator with
#polynomial coefficients, in terms of z and the
#differentiation w.r.t. z, Dz, that annihilates
#the (ordinary) generating function
#of the sequence
#"number of words in {1, ...m} of length n that has
#as many factors (substrings) that are w1 and
#factors that are w2 equals 
#such that the constant term of s is the generating
#function for words in {1,..,m}, according to length that
#have the same number of occurences of the factor w1
#and w2.
#For example, to completely solve Richard Stanley's
#Amer. Math. Monthly problem #11610 (page 937
# of v. 118 (No. 10), [Dec. 2011]), please type:
#DOaz(2,[1,1],[1,2],t,Dt);
DOaz:=proc(m,w1,w2,z,Dz) local mu,gu,s:
option remember:
if not (type(w1,list) and type(w2,list) and 
{op(w1)} minus {seq(i1,i1=1..m)}={} and
{op(w2)} minus {seq(i1,i1=1..m)}={} and
w1<>w2) then
 print(`bad input`):
 RETURN(FAIL):
fi:

mu:=RSct(m,w1,w2,z,s)/s:
gu:=AZc(mu,s,z,Dz):
if gu<>0 then
 RETURN(gu[1]):
else
 RETURN(FAIL):
fi:
end:


#ProveGF(m,w1,w2,z,F): Given a positive integer
#m and two words (lists),w1,w2, of the same length in the
#alphabet {1, ...,m}, and an expression F that
#depends on z,
#proves (rigorously!), using
#the continuous Almkvist-Zeilberger algorithm,
#that F is the ordinary generating function (using 
#the variable z) of the sequence
#"number of words in {1, ...m} of length n that has
#as many factors (substrings) that are w1 and
#factors that are w2 equals 
#For example, to completely solve Richard Stanley's
#Amer. Math. Monthly problem #11610 (page 937
# of v. 118 (No. 10), [Dec. 2011]), please type:
#ProveGF(2,[1,1],[1,2],t,1/(1-t)/2+(1+2*t)/sqrt((1-t)*(1-2*t)*(1+t+2*t^2))/2);
ProveGF:=proc(m,w1,w2,z,F) local mu,gu,s,richard,Dz,i,mu1,F1,Deg1,i1,
ku1,ku2:
mu:=RSct(m,w1,w2,z,s):
gu:=AZc(mu/s,s,z,Dz):
if gu=0 then
 RETURN(FAIL):
fi:

gu:=gu[1]:
Deg1:=degree(gu,Dz)+degree(gu,z):
F1:=taylor(F,z=0,Deg1+1):

mu1:=taylor(mu,z=0,Deg1+1):
ku1:=[seq(coeff(expand(coeff(mu1,z,i1)),s,0),i1=0..Deg1)]:
ku2:=[seq(coeff(F1,z,i1),i1=0..Deg1)]:

if ku1<>ku2 then
print(`Boundary conditions do not match`):
print(ku1,ku2):
 RETURN(false):
fi:

richard:=add(coeff(gu,Dz,i)*diff(F,z$i),i=1..degree(gu,Dz))+
coeff(gu,Dz,0)*F:
richard:=normal(simplify(richard)):
if richard<>0 then
 print(`Diff. eq. not satisfied`):
 RETURN(false):
fi:
true:
end:





#RECaz(m,w1,w2,n,N): Given a positive integer
#m and two words (lists),w1,w2, of the same length in the
#alphabet {1, ...,m}, 
#uses the continuous Almkvist-Zeilberger algorithm,
#to find a linear recurrence operator with
#polynomial coefficients, ope(n,N) in terms of n and the
#shift operator in n, N, that annihilates
#the sequence
#"number of words in {1, ...m} of length n that has
#as many factors (substrings) that are w1 and
#factors that are w2 equals 
#such that the constant term of s is the generating
#function for words in {1,..,m}, according to length that
#have the same number of occurences of the factor w1
#and w2. It also gives the initial values.
#For example, to get a recurrence operator for
#Richard Stanley's
#Amer. Math. Monthly problem #11610 (page 937
# of v. 118 (No. 10), [Dec. 2011]), please type:
#RECaz(2,[1,1],[1,2],n,N);
RECaz:=proc(m,w1,w2,n,N) local ope,z,Dz,ORDER,lu,i1:
option remember:
ope:=DOaz(m,w1,w2,z,Dz):
ope:=DiffToRec1(ope,z,Dz,n,N):
ORDER:=degree(ope,N):
ope:=add(factor(coeff(ope,N,i1))*N^i1,i1=0..ORDER):
lu:=RSseq(m,w1,w2,ORDER):
[ope,lu]:
end:



#RECg(m,w1,w2,n,N,MaxC): Given a positive integer
#m and two words (lists),w1,w2, of the same length in the
#alphabet {1, ...,m}, 
#uses guessing 
#to find a linear recurrence operator with
#polynomial coefficients, ope(n,N) in terms of n and the
#shift operator in n, N, that annihilates
#the sequence
#"number of words in {1, ...m} of length n that has
#as many factors (substrings) that are w1 and
#factors that are w2 equals 
#such that the constant term of s is the generating
#function for words in {1,..,m}, according to length that
#have the same number of occurences of the factor w1
#and w2. It also gives the initial values.
#For example, to get a recurrence operator for
#Richard Stanley's
#Amer. Math. Monthly problem #11610 (page 937
# of v. 118 (No. 10), [Dec. 2011]), please type:
#RECg(2,[1,1],[1,2],n,N,10);
RECg:=proc(m,w1,w2,n,N,MaxC) local gu,K,ope:
option remember:
K:=max(trunc(MaxC^2/4+12),7*MaxC+6):

gu:=RSseq(m,w1,w2,K):
ope:=Findrec(gu,n,N,MaxC):

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

[ope,[op(1..degree(ope,N),gu)]]:
end:

#Khalek(ope1,ope2,n,N): given linear recurrence operators ope1
#and (monin in N) ope2 in n and N (the forward shift operator in N),
#outputs operators Q(N,n) and R(N,n) such that
#ope1(N)=Q(N,n)ope2(n,N)+R(N,n) where the order of Q(N,n)(i.e. degree in N)
#is the difference in the oders of ope1 and ope2 and the order of
#R(N,n) is less than the order of ope2
#For example, try
#Khalek(n*N^2+n/(n+1)*N,N+1/n,n,N);
Khalek:=proc(ope1,ope2,n,N) local gu,lu,ORDER1,ORDER2,ope1A,i1:
ORDER1:=degree(ope1,N):
ORDER2:=degree(ope2,N):

if coeff(ope2,N,ORDER2)<>1 then
 RETURN(FAIL):
fi:

if degree(ope1,N)<degree(ope2,N) then
 RETURN([0,ope1]):
fi:

lu:=coeff(ope1,N,ORDER1)*N^(ORDER1-ORDER2):
ope1A:=expand(ope1-lu*subs(n=n+ORDER1-ORDER2,ope2)):
ope1A:=add(normal(coeff(ope1A,N,i1))*N^i1,i1=0..degree(ope1A,N)):
gu:=Khalek(ope1A,ope2,n,N):
[lu+gu[1],gu[2]]:
end:







#ALGg(m,w1,w2,K,x,P): Given a positive integer
#m and two words (lists),w1,w2, of the same length in the
#alphabet {1, ...,m},w1, w2, and a pos. integer K
#uses guessing with K data points to find an algebraic equation
#in x and P
#satisfied by the ordinary generating function, in x,
#let's call it P,
#of the sequence
#"number of words in {1, ...m} of length n that has
#as many factors (substrings) that are w1 and
#factors that are w2 equals 
#such that the constant term of s is the generating
#function for words in {1,..,m}, according to length that
#have the same number of occurences of the factor w1
#and w2. 
#For example, to get an algebraic equation for
#Richard Stanley's
#Amer. Math. Monthly problem #11610 (page 937
# of v. 118 (No. 10), [Dec. 2011]), please type:
#ALGg(2,[1,1],[1,2],40,x,P);
ALGg:=proc(m,w1,w2,K,x,P) local gu,i1:
option remember:

if not type(m,integer) and m>1 and type(K,integer) and K>0 and
type(x,symbolic) and type(P,symbolic) then
  print(`bad input`):
 RETURN(FAIL):
fi:
 
if not (type(w1,list) and type(w2,list) and 
{op(w1)} minus {seq(i1,i1=1..m)}={} and
{op(w2)} minus {seq(i1,i1=1..m)}={} and
w1<>w2) then
 print(`bad input`):
 RETURN(FAIL):
fi:


gu:=RSseq(m,w1,w2,K):
gu:=[1,op(gu)]:
Empir(gu,x,P):
end:

rev11:=proc(w) local i:
[seq(w[nops(w)-i+1],i=1..nops(w))]:
end:

rev1:=proc(S) local w:
{seq(rev11(w),w in S)}:
end:

#Images1(m,S): Given a set of words S in the alphabet
#{1, ...,m}, finds the set of images under permuting
#the letters and reversing the words
#For example, try:
#Images1(2,{[1,1],[1,2]});
Images1:=proc(m,S) local mu,gu1,pi,gu11,i,gu2:
mu:=permute(m):
gu1:={seq(subs({seq(i=pi[i],i=1..m)},S),pi in mu)}:
gu2:={seq(rev1(gu11),gu11 in gu1)}:
gu1 union gu2:
end:


#Zugot(m,k): Representatives of sets of pairs of words of
#length k in the alphabet {1, ..., m}. For example, try:
#Zugot(2,2);
Zugot:=proc(m,k) local gu,mu:

mu:=Milim(m,k):
mu:=choose(mu,2);

gu:={}:

while mu<>{} do
gu:=gu union {mu[1]}:
mu:=mu minus Images1(m,mu[1]):
od:
gu:
end:








#Info1(m,w1,w2,n,N,P,t,K): Given a pos. integer m,
#and two words w1 and w2 of the same length in the
#alphabet {1, ..., m}, and symbols n and N
#P and t.  It also inputs a positive integer
#K (the number of terms for guessing the alg. relation, see below)
#We are interested in the integer sequence
#a(n):=number of n-lettered words in the alphabet {1, ...,m}
#that has as many substrings (alias factors, 
#alias consecutive subwords) that are w1 as there are w2.
#The output this procedure  is a list consisting
#(i) a (conjectured) recurrence operator ope(n,N) annihilating
#the sequence a(n) together with the initial values 
#(starting at n=1) providing
#a full holonomic (P-recursive) description of the sequence
#(ii) The minimal polynomial F(P,x) such that
#F(P(x),x)=0 where P(x) is the ordinary generating function of
#the sequence: i.e.,
#P(x)=Sum(a[n]*x^n,n=0..infinity);
#(iii) the "explicit" expression for P(x) in terms of radicals.
#For example, for Richard Stanley's Amer. Math. Monthly Problem
#11610 (118(10)[Dec. 2011], p. 937)
#type
#Info1(2,[1,1],[1,2],n,N,P,t,40);
Info1:=proc(m,w1,w2,n,N,P,t,K) local MaxC,opeD,REC1,ku:
option remember:
opeD:=RECaz(m,w1,w2,n,N):
MaxC:=degree(opeD[1],N)+degree(numer(opeD[1]),n):

REC1:=RECg(m,w1,w2,n,N,MaxC):

ku:=ALGg(m,w1,w2,K,t,P):

if ku=0 then

print(`K =`, K, `is too small, make it bigger`):
 RETURN(FAIL):
fi:

 [[m,{w1,w2}],REC1,ku]:

end:


#Info2(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m,
#and two words w1 and w2 of the same length in the
#alphabet {1, ..., m}, and symbols n and N
#P and t.  It also inputs a positive integer
#K (the number of terms for guessing the alg. relation, see below)
#We are interested in the integer sequence
#a(n):=number of n-lettered words in the alphabet {1, ...,m}
#that has as many substrings (alias factors, 
#alias consecutive subwords) that are w1 as there are w2.
#The output this procedure  is a list consisting
#(i) a (conjectured) recurrence operator ope(n,N) annihilating
#the sequence a(n) together with the initial values 
#(starting at n=1) providing
#a full holonomic (P-recursive) description of the sequence
#(ii) The minimal polynomial F(P,x) such that
#F(P(x),x)=0 where P(x) is the ordinary generating function of
#the sequence: i.e.,
#P(x)=Sum(a[n]*x^n,n=0..infinity);
#(iii) The asymptotics of a(n) to order M
#For example, for Richard Stanley's Amer. Math. Monthly Problem
#11610 (118(10)[Dec. 2011], p. 937)
#type
#Info2(2,[1,1],[1,2],n,N,P,t,40,3);
Info2:=proc(m,w1,w2,n,N,P,t,K,M) local MaxC,opeD,REC1,ku,lu,C,vu,i1,C1:
option remember:
opeD:=RECaz(m,w1,w2,n,N):
MaxC:=degree(opeD[1],N)+degree(numer(opeD[1]),n):

REC1:=RECg(m,w1,w2,n,N,MaxC):

ku:=ALGg(m,w1,w2,K,t,P):

if ku=0 then

print(`K =`, K, `is too small, make it bigger`):
 RETURN(FAIL):
fi:

lu:=AsyC(REC1[1],n,N,M,REC1[2],2000);

if lu<>FAIL then
C:=lu[1]:
if lu[3]>8 then
C1:=identify(lu[1]):
else
 C1:=0:
fi:
lu:=C*lu[2]:
fi:

vu:=SeqFromRec(REC1[1],n,N,evalf(REC1[2]),2000):
vu:=evalf([seq(vu[i1]/subs(n=i1,lu),i1=2000-10..2000)]):

 [[m,{w1,w2}],REC1,ku,lu,vu,C1,C]:
end:


#TerseBook(m,k,n,N,P,t,K): inputs positive integers m and k
#and outputs a set of lists consiting of 
#the input and output of Info1(m,k,n,N,P,t,K,M)
#(qv) for all members of Zugot(m,k) (qv).
#For example, try:
#TerseBook(2,2,n,N,P,t,40);
TerseBook:=proc(m,k,n,N,P,t,K) local gu,i:
option remember:
gu:=Zugot(m,k):
{seq(Info2(m,op(gu[i]),n,N,P,t,K),i=1..nops(gu))}:
end:



#HSP(ope,n,N): the highest singularity of the recurrence
#operator ope(n,N). For example, try:
#HSP((n-3)*(n-4)*N+1,n,N);
HSP:=proc(ope,n,N) local ope1,pol,gu,i,gu1:
ope1:=numer(normal(ope)):
pol:=lcoeff(ope1,N):
gu:=[solve(pol,n)]:
gu1:={}:
for i from 1 to nops(gu) do
 if type(gu[i],integer) and gu[i]>=0 then
   gu1:=gu1 union {gu[i]}:
 fi:
od:

max(op(gu1)):
end:






#Mishpat(m,w1,w2,n,N,P,t,K,M): Given a pos. integer m,
#and two words w1 and w2 of the same length in the
#alphabet {1, ..., m}, and symbols n and N
#P and t.  It also inputs a positive integer
#K (the number of terms for guessing the alg. relation, see below)
#We are interested in the integer sequence
#a(n):=number of n-lettered words in the alphabet {1, ...,m}
#that has as many substrings (alias factors, 
#alias consecutive subwords) that are w1 as there are w2.
#The output this procedure  is a verbose theorem and sketch of proof
#about the sequence a(n)
#For example, for Richard Stanley's Amer. Math. Monthly Problem
#11610 (118(10)[Dec. 2011], p. 937)
#type
#Mishpat(2,[1,1],[1,2],n,N,P,t,40,3);
Mishpat:=proc(m,w1,w2,n,N,P,t,K,M) local gu,i,mu,ope,i1,a,d1:

gu:=Info2(m,w1,w2,n,N,P,t,K,M):

print(`Proposition: Let a(n) be the number of n-lettered words`):
print(`in the alphabet`, {seq(i,i=1..m)}):
print(`With as many occurrences of the substring (consecutive subword)`, w1):
print(`as those of`, w2):
print(``):
print(`Let P(t), be the ordinary generating function of a(n)`):
print(P(t)=Sum(a(n)*t^n,n=0..infinity)):

mu:=gu[3]:

if degree(mu,P)=2 then

print(`Then P=P(t) satisfies the quadratic equation`):
print(mu=0):
else
print(`Then P=P(t) satisfies the algebaric equation`):
print(mu=0):
fi:

ope:=gu[2]:


print(`Not only that, a(n)  safisfies the following linear recurrence`):
d1:=degree(ope[1],N):
print(a(n+degree(ope[1],N))=
-add(coeff(ope[1],N,d1-i1)*a(n+d1-i1),i1=1..d1)):
print(``):
print(`subject to the initial conditions`):
print(seq(a[i1]=ope[2][i1],i1=1..nops(ope[2]))):

if gu[4]<>FAIL then
 print(`Finally, the asymptotics of a(n) to order`,M, `is : `):
 print(gu[4]):
 print(`Note that everything is rigorous except the constant in front`,gu[7]):
 print(`that is a mere non-rigorous estimate`):
 if gu[6]<>0 and gu[6]<>gu[7] then
  print(` Maple conjectures that its exact value`):
  print(`is equal to`, gu[6]):
 fi:

fi:

print(``):

print(`Sketch of proof: The above has been proved completely rigorously`):
print(`internally, and I will spare you the details. If you insist you`): 
print(` can get them yourself by tracing the code, or else donate `):
print(` 50 dollars  to the Society for the Liberation of Computerkind.`):
print(`QED. `):
print(``):
print(`For the sake of people who want to use these formulas here`):
print(`are the quardaric equation, recurrence operator, and asymptotics`):
print(`in Maple input format`):
lprint(mu,ope,gu[4]):

print(``):
print(`For the sake of Sloane, here are the first 50 terms of the sequece`):
print([1,op(RSseq(m,w1,w2,50))]):

end:


#Mishpat1(m,w1,w2,n,N,P,t,K,M,k): 
#Like Mishpat(m,w1,w2,n,N,P,t,K,M,k) but with the Proposition
#numbered k
#For example, for Richard Stanley's Amer. Math. Monthly Problem
#11610 (118(10)[Dec. 2011], p. 937)
#type
#Mishpat1(2,[1,1],[1,2],n,N,P,t,40,3,1);
Mishpat1:=proc(m,w1,w2,n,N,P,t,K,M,k) local gu,i,mu,ope,i1,a,d1:

gu:=Info2(m,w1,w2,n,N,P,t,K,M):

print(`Proposition Number `, k, `: Let a(n) be the number of n-lettered words`):
print(`in the alphabet`, {seq(i,i=1..m)}):
print(`With as many occurrences of the substring (consecutive subword)`, w1):
print(`as those of`, w2):
print(``):
print(`Let P(t), be the ordinary generating function of a(n)`):
print(P(t)=Sum(a(n)*t^n,n=0..infinity)):

mu:=gu[3]:

if degree(mu,P)=2 then

print(`Then P=P(t) satisfies the quadratic equation`):
print(mu=0):
else
print(`Then P=P(t) satisfies the algebaric equation`):
print(mu=0):
fi:

ope:=gu[2]:


print(`Not only that, a(n)  safisfies the following linear recurrence`):
d1:=degree(ope[1],N):
print(a(n+degree(ope[1],N))=
-add(coeff(ope[1],N,d1-i1)*a(n+d1-i1),i1=1..d1)):
print(``):
print(`subject to the initial conditions`):
print(seq(a[i1]=ope[2][i1],i1=1..nops(ope[2]))):

if gu[4]<>FAIL then
 print(`Finally, the asymptotics of a(n) to order`,M, `is : `):
 print(gu[4]):
 print(`Note that everything is rigorous except the constant in front`,gu[7]):
 print(`that is a mere non-rigorous estimate`):
 if gu[6]<>0 and gu[6]<>gu[7] then
  print(` Maple conjectures that its exact value`):
  print(`is equal to`, gu[6]):
 fi:

fi:

print(``):

print(`Sketch of proof: The above has been proved completely rigorously`):
print(`internally, and I will spare you the details. If you insist you`): 
print(` can get them yourself by tracing the code, or else donate `):
print(` 50 dollars  to the Society for the Liberation of Computerkind.`):
print(`QED. `):
print(``):
print(`For the sake of people who want to use these formulas here`):
print(`are the quardaric equation, recurrence operator, and asymptotics`):
print(`in Maple input format`):
lprint(mu,ope,gu[4]):

print(``):
print(`For the sake of Sloane, here are the first 50 terms of the sequece`):
print([1,op(RSseq(m,w1,w2,50))]):

end:







#VerboseBook(m,k,n,P,t,K,M): inputs positive integers m and k
#and outputs a web-book with generating functions
#and recurrences for sequences enumerating
#analogous quantities to Richard Stanley's Monthly
#Problem 11610 (118(10)[Dec. 2011], p. 937)
#concerning all pairs of words in {1, ...,m} 
#of length k (up to equivalence). 
#m is the size of the alphabet, k is the length of the
#distinguished words, n,P,t is the symbol for displaying
#K and M are as in Info2 
#For example try:
#VerboseBook(2,2,n,P,t,40,4);
VerboseBook:=proc(m,k,n,P,t,K,M) local gu,i1,t0,N:
gu:=Zugot(m,k):
t0:=time():
print(`Analogs of the Richard Stanley Amer. Math. Monthly Problem  11610`):
print(`for ALL pairs of words of length`, k, `in an alphabet of`,m):
print(`letters.`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for i1 from 1 to nops(gu) do
 Mishpat1(m,op(gu[i1]),n,N,P,t,K,M,i1):
od:

print(`It took`, time()-t0, `seconds to generate this book. `):

end:




#CheckRec(L,ope,n,N): Given a sequence of numbers L and a
#recurrence operator ope(n,N) checks whether indeed
#L is annihilated by ope(n,N). For example, try:
#CheckRec([1,2,4,8,16],N-2,n,N);
CheckRec:=proc(L,ope,n,N) local n1,i1:
evalb({0}=
{seq(
add(subs(n=n1,coeff(ope,N,i1))*L[n1+i1],i1=0..degree(ope,N)),
n1=1..nops(L)-degree(ope,N))}):
end:


#CheckAsy(a,n,L,eps): checks whether the asympotic formula  a, an
#expression in n, is a good one for the sequence L. With tolerance
#eps. For example, try:
#CheckAsy(2^n,n,[seq(2^i,i=1..100)],1/100);
CheckAsy:=proc(a,n,L,eps) local n1,gu:

gu:=
min(
seq(evalf(abs(L[n1]/subs(n=n1,a)-1.)),n1=trunc(nops(L)/2)..nops(L))):

if gu<eps then
RETURN(true):
else
RETURN(false):
fi:

end: