######################################################################
##FunEq: Save this file as   FunEq                                   #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read    FunEq<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  Dec. 25, 2014`):
print(`This version:  Jan. 5, 2015`):
print(` This is FunEq  `):
print(`A Maple package that uses an empirical-yet-fully-rigorous method to solve functional equations of the form`):
print(``):
print(` F(P(x,t), P(x,1),x,t))=0 . `):
print(``):

print(`It accompanyies the article (or rather the article accompanyies this Maple package!, the tail is wagging the dog!)`):
print(``):
print(` An Empirical Method for Solving (rigorously!) Algebraic Functional Equations Of the Form`):
print(` F(P(x,t), P(x,1),x,t))=0 `):
print(`by Ira M. Gessel and Doron Zeilberger`):
print(`published exclusively in the Personal Journal of S.B. Ekhad and D. Zeilberger: http://www.math.rutgers.edu/~zeilberg and arxiv.org`):
print(`and Ira Gessel's website, but NOT in a "regular" journal. `):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

print(`---------------------------------------`):
 print(`For a list of the Story procedures, that output articles or books, type ezraV();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):


print(`---------------------------------------`):
 print(`For a list of the SCHUTZENBERGER procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


with(combinat):



ezra01:=proc()

if args=NULL then
 print(` The  procedures dealing with only Beta(0,1) trees are: `):
   print(` ,  Bab, BetaEq01, CheckFunEq01, FunEq01, PreBetaEq,, PreT01, T01 `):

else
ezra(args):
fi:

end:

ezraV:=proc()

if args=NULL then
 print(` The  Verbose procedures are: `):
 print(` InfoV, PrintFE, SeferAlpha, SeferBeta, SeferJW `):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` AsyAnal,  empir, Empir, FEtoAEslow, OneStep, OneStepBab, OneStepBabK, OneStepK,  Zinn `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: AabFE, AEtoSeq,  BabFE, FEtoAE, FEtoAE1, FEtoAE1wp, FEtoRec1, FirstK, GenAabFE, GenBabFE, GenjwFE, Info, jwFE, MBM   `):
 print(` `):

elif nops([args])=1 and op(1,[args])=AEtoSeq then
print(`AEtoSeq(F,P,x,a0,K): given a polynomial F in the variables  P and x, where a0 is the constant term`):
print(`finds the list of coefficients of x^i from i=1, to K of the unique formal power series P(x)`):
print(` satisfying the algebraic equation `):
print(` F(P(x),x)=0 and P(0)=a0. Try: `):
print(` AEtoSeq(C-1-x*C^2,C,x,1,20) ; `):

elif nops([args])=1 and op(1,[args])=AsyAnal then
print(`AsyAnal(L,N,n): given a list of numbers L, of length k, say, whose i-th entry is f(N+i-1) for i from 1 to nops(L), conjectures`):
print(`an asymptotic approximate expression for f(m) in the form a[0]+a[1]/n+...+ a[k]/n^(k-1). Try:`):
print(`AsyAnal([seq(2+6/i+11/i^2,i=98..100)],98,n);`):



elif nops([args])=1 and op(1,[args])=Bab then
print(`Bab(a,b,y,K): The first K terms in the sequence of weight-enumetators, according to the weight y^(label of root)`):
print(`of Beta(a,b) description trees, that satisies the functional equation`):
print(`F=x+x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F)`):
print(`Following the article:`):
print(` Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer.  `):
print(`Description trees for some families of planar maps." Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.`):
print(` http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz `):
print(` Try: `):
print(`Bab(0,1,y,10);`):

elif nops([args])=1 and op(1,[args])=AabFE then
print(`AabFE(a,b,x,y,fxy,fx1): the functional equation for A(a,b) description trees, given`):
print(`as a function of x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1)), `):
print(` the functional equation itself is `):
print(` f(x,y)=F(x,y,f(x,y),f(x,1)); Try: `):
print(`AabFE(0,1,x,y,fxy,fx1);`):
elif nops([args])=1 and op(1,[args])=BabFE then
print(`BabFE(a,b,x,y,fxy,fx1): the functional equation for B(a,b) description trees, given`):
print(`as a function of x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1)), `):
print(` the functional equation itself is `):
print(` f(x,y)=F(x,y,f(x,y),f(x,1)); Try: `):
print(`BabFE(0,1,x,y,fxy,fx1);`):

elif nops([args])=1 and op(1,[args])=BetaEq01 then
print(` BetaEq01(P,x,z,K) : Guessing an algebraic equation satisfied by the generating function of Pre-Beta-01-trees, using K terms, try: `):
print(` BetaEq01(P,x,z,30); `):

elif nops([args])=1 and op(1,[args])=CheckFunEq01 then
print(`CheckFunEq01(F,x,z,K): Checks that F satisfies  the Functional Equation for F(x,z) for Pre-01-trees up to K terms, Try:`):
print(` gu:=PreT01(20,z): gu:=add(gu[i]*x^i,i=1..20):  CheckFunEq01(gu,x,z,20); `):

elif nops([args])=1 and op(1,[args])=empir then
print(`empir(gu,degx,degP,x,P) empirically finds an algebraic equation`):
print(` F(P(x),x)=0 of degree degP in P(x) and degx in x for `):
print(`formal power series P(x):=sum_i gu[i]*x^i, where gu is a list`):
print(`for example empir([seq(1,i=1..20)],1,1,x,P) should yield`):
print(` P-xP-1=0`):
 
elif nops([args])=1 and op(1,[args])=Empir then
print(`Empir(gu,x,P) empirically finds an algebraic equation`):
print(` F(P(x),x)=0 for the  `):
print(`formal power series P(x):=sum_{i=0} gu[i]*x^i, where gu is a list`):
print(`for example Empir([seq(1,i=1..20)],x,P) should yield`):
print(` P-xP-1=0`):
print(`If there is not enough data, it reurns 0. You should then`):
print(`try it again with a longer sequence`):
 

elif nops([args])=1 and op(1,[args])=FEtoAE then
print(`FEtoAE(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))`):
print(`guesses, using the first K terms, an algebraic equation for f(x,y), denoted by fxy. Try:`):
print(`FEtoAE(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);`):

elif nops([args])=1 and op(1,[args])=FEtoAEslow then
print(`STUPID SLOW WAY OF DOING FEtoAE(F,x,y,fxy,fx1,K), ONLY FOR CHECKING PURPOSES. `):
print(`FEtoAEslow(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))`):
print(`guesses, using the first K terms, an algebraic equation for f(x,y), denoted by fxy. Try:`):
print(`FEtoAEslow(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);`):



elif nops([args])=1 and op(1,[args])=FEtoAE1 then
print(`FEtoAE1(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))`):
print(`guesses, using the first K terms, an algebraic equation for f(x,1).`):
print(`If insufficient data, it returns FAIL. Try:`):
print(`FEtoAE1(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);`):


elif nops([args])=1 and op(1,[args])=FEtoAE1wp then
print(`FEtoAE1wp(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))`):
print(`guesses, using the first K terms, an algebraic equation for f(x,1). `):
print(`IT THEN GOES ON TO PROVE IT!`):
print(`If it fails it tells you so. Try:`):
print(`FEtoAE1wp(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);`):

elif nops([args])=1 and op(1,[args])=FEtoAE1wpV then
print(`FEtoAE1wpV(F,x,y,fxy,fx1,K,f):  Verbose form of FEtoAE1wp(F,x,y,fxy,fx1,K) (q.v.) `):
print(`Try:`):
print(`FEtoAE1wpV(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30,f);`):

elif nops([args])=1 and op(1,[args])=FEtoRec1 then
print(`FEtoRec1(F,x,y,fxy,fx1,n,N,C): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))`):
print(`guesses a recurrence operator in terms of n and the shift operator N of comlexity <=C annihilating the`):
print(`coefficents of f(x,1). Try: `):
print(`FEtoRec1(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,n,N,6); `):

elif nops([args])=1 and op(1,[args])=FirstK then
print(`FirstK(F,x,y,fxy,fx1,K): The list of coeficients of x^i, from i=1 to i=K `):
print(`of the formal power series in x,y, that satisfies the functional equation`):
pritn(` f(x,y)=F(x,y,f(x,y),f(x,1)). Try: `):
print(` FirstK(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,10); `):

elif nops([args])=1 and op(1,[args])=FunEq01 then
print(`FunEq01(F,x,z): The Functional Equation for F(x,z): the genrating function for Pre-01-trees where the weight is x^(#vertices)*z^(label of root), Try:`):
print(`FunEq01(F,x,z)`):

elif nops([args])=1 and op(1,[args])=GenAabFE then
print(`GenAabFE(a,b,A,B,x,y,fxy,fx1), a generalization of AabFE(a,b,x,y,fxy,fx1) with parameters A,B`):
print(`GenAabFE(a,b,1,1,x,y,fxy,fx1) is the same as AabFE(a,b,xy,fxy,fx1): Try:`):
print(`GenAabFE(1,1,1,1,x,y,fxy,fx1);`):

elif nops([args])=1 and op(1,[args])=GenBabFE then
print(`GenBabFE(a,b,A,B,C,x,y,fxy,fx1), a generalization of BabFE(a,b,x,y,fxy,fx1) with parameters A,B,C`):
print(`GenBabFE(a,b,1,1,1,x,y,fxy,fx1) is the same as BabFE(a,b,xy,fxy,fx1): Try:`):
print(`GenBabFE(0,1,1,1,1,x,y,fxy,fx1);`):

elif nops([args])=1 and op(1,[args])=GenjwFE then
print(`GenjwFE(x,y,fxy,fx1,A,B): The Generalization of the Functional Equation for Two-Stack-Sortable Permutations`):
print(`with two extra parameters, A,B. When A=1, B=1 it is the original. Try:`):
print(`GenjwFE(x,y,fxy,fx1,2,3);`):

elif nops([args])=1 and op(1,[args])=Info then
print(`Info(F,x,y,fxy,fx1,K,n,N,MaxC,L0,L1):  `):
print(` Inputs: `):
print(` (i) A polynomial F in the variables x,y,fxy,fx1, representing a functional equation F(f(x,y),f(x,1),x,y)=0 `):
print(` (ii) variables x,y,fx1,fx1, for expressing the above-mentioned Functional Equation `):
print(`(iii) a positive integer K indicating how many terms to use for guessing the algebraic equation`):
print(`(iv) symbols n and N for guessing a linear recurrence operator ope(n,N) annihilating the sequence`):
print(`(v): MaxC: a positive integer MaxC indicating the maximum complexity of the recurrence`):
print(`(vi) positive integers L0 and L1, `):
print(` outputs (in the best possible case) a list [eq,ope,ListOfFirstL0 terms, L1-th term] `):
print(` Try: `):
print(` Info(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,20,n,N,10,30,1000); `):

elif nops([args])=1 and op(1,[args])=InfoV then
print(`InfoV(F,x,y,fxy,fx1,K,n,N,MaxC,L0,L1):  verbose form of Info(F,x,y,fxy,fx1,K,n,N,MaxC,L0,L1) (q.v.)`):
print(` Try: `):
print(` InfoV(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,20,n,N,10,30,100); `):

elif nops([args])=1 and op(1,[args])=jwFE then
print(`jwFE(x,y,fxy,fx1): The Functional Equation for Two-Stack-Sortable Permutations`):

elif nops([args])=1 and op(1,[args])=MBM then
print(`MBM(F,x,y,fxy,fx1): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))`):
print(`finds an algebraic equation for f(x,1).`):
print(`using Mireille's Bousquet-Melou's amazing double disriminant. Try:`):
print(`MBM(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);`):

elif nops([args])=1 and op(1,[args])=OneStep then
print(`OneStep(F,x,z): given a formal power series F in x,z, performs one iteration in`):
print(`F->x+x*T(F), where T(f)=(subs(z=1,f)-z^2*f)/(1-z)`):
print(`Try:`):
print(`OneStep(x,x,z);`):


elif nops([args])=1 and op(1,[args])=OneStepBab then
print(`OneStepBab(F,x,y,a,b): given a formal power series F in x,y, performs one iteration in`):
print(`the evolution of Beta(a,b)-Description trees according to the article:`):
print(`F->x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F)`):
print(`Following the article:`):
print(` Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer.  `):
print(`Description trees for some families of planar maps." Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.`):
print(` http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz `):
print(` Try: `):
print(`OneStepBab(x,x,y,0,1);`):

elif nops([args])=1 and op(1,[args])=OneStepBabK then
print(`OneStepBabK(F,x,y,a,b,K): given a formal power series F in x,y, performs one iteration in`):
print(`the evolution of Beta(a,b)-Description trees according to the article:`):
print(`F->x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F)`):
print(`and truncates it to order K in x. `):
print(`Following the article:`):
print(` Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer.  `):
print(`Description trees for some families of planar maps." Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.`):
print(` http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz `):
print(` Try: `):
print(`OneStepBabK(x,x,y,0,1,10);`):

elif nops([args])=1 and op(1,[args])=OneStepK then
print(`OneStepK(F,x,z): given a formal power series F in x,z, performs one iteration in`):
print(`F->x+x*T(F), where T(f)=(subs(z=1,f)-z^2*f)/(1-z), and truncates it to order K in x`):
print(`Try:`):
print(`OneStepK(x,x,z,10);`):

elif nops([args])=1 and op(1,[args])=PreBetaEq01 then
print(` PreBetaEq01(P,x,z,K) : Guessing an algebraic equation satisfied by the generating function of Pre-Beta-01-trees, using K terms, try: `):
print(` PreBetaEq01(P,x,z,30); `):


elif nops([args])=1 and op(1,[args])=PreT01 then
print(`PreT01(K,z): The first K terms in the sequence of the generating function, according to the label of the root, of Pre-Beta01 trees`):
print(` Try: `):
print(` PreT01(10,z); `):

elif nops([args])=1 and op(1,[args])=PrintFE then
print(`PrintFE(F,x,y,fxy,fx1,f): spells hour the functional equation in human form, try:`):
print(`PrintFE(jwFE(x,y,fxy,fx1),x,y,fxy,fx1,f); `):

elif nops([args])=1 and op(1,[args])=SeferAlpha then
print(`SeferAlpha(K,C,L0,L1,g): inputs positive integers K,C,L0,L1 and g and outputs an article about`):
print(` the enumeration of  Alpha(a,b) Description trees for 0<=a,b<=g. For each it also `):
print(` outputs the linear recurrence equation with polynomial coefficients satisfied by the enumerating sequence `):
print(` if its ORDER+DEGREE<=C. `):
print(` It also outputs the set of all  successful `):
print(` pairs [a,b] for which K data points suffice. For example, Try: `):
print(` SeferAlpha(20,10,30,100,2); `):

elif nops([args])=1 and op(1,[args])=SeferBeta then
print(`SeferBeta(K,C,L0,L1,g): inputs positive integers K,C,L0,L1 and g and outputs an article about`):
print(` the enumeration of  Beta(a,b) Description trees for 0<=a,b<=g. For each it also `):
print(` outputs the linear recurrence equation with polynomial coefficients satisfied by the enumerating sequence `):
print(` if its ORDER+DEGREE<=C. `):
print(` It also outputs the set of all  successful `):
print(` pairs [a,b] for which K data points suffice. For example, Try: `):
print(` SeferBeta(20,10,30,100,2); `):

elif nops([args])=1 and op(1,[args])=SeferJW then
print(`SeferJW(K,C,L0,L1): inputs positive integers K,C,L0,L1 and g and outputs an article about`):
print(`the enumeration of  2-stack-sortable permutations`):
print(`Try: `):
print(`SeferJW(20,10,30,100);`):

elif nops([args])=1 and op(1,[args])=T01 then
print(`T01(K,z): The first K terms in the sequence of the generating function, according to the label of the root, of Beta01 trees`):
print(` Try: `):
print(` T01(10,z); `):

elif nops([args])=1 and op(1,[args])=Zinn then 
print(`Zinn(resh): Zinn-Justin's method to estimate`):
print(`the C,mu, and theta such that`):
print(`resh[i] is appx. Const*mu^i*i^theta`):
print(`For example, try:`):
print(`Zinn([seq(5*i*2^i,i=1..30)]);`):

print(``):


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


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

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 SCHUTZENBERGER
ezraS:=proc()
if args=NULL then
 print(`Contains the following procedures:`):
 print(`difftorec, algtodiff,genseq,genseqnew,empir,Findrec,mfindrec,logalg`):
 print(`algprod,newts,newt,pol_in_terms_of_p,decomp,decompx,findrec2`):
 print(` Empir , qfac1,qbin1, qfindrec `):
 print(`a specific procedure, type ezra(procedure_name)`):
print(`Warning: q is a global variable, do not use it!`):
fi:
 
 
if nops([args])=1 and op(1,[args])=qfac1 then
print(`qfac1(n): (1-q)*...*(1-q^n) `):
fi:
 
if nops([args])=1 and op(1,[args])=qbin1 then
print(`qbin1(n,k): qfac1(n)/qfac1(k1)/qfac1(n-k) `):
fi:
 
if nops([args])=1 and op(1,[args])=`decompx` then
 
print(`decompx(pol,N,x): given a polynomial in the variables N and x`):
print(`investigates whether, when viewed as a polynomial in N`):
print(` it can be written as the alg. product (i.e. its roots are`):
print(`the direct product of the sets of the roots of its " factors")`):
print(`of a quadratic in N`):
print(`of the form N^2+2*a*x*N-1 where a is an (algebraic) pure number`):
print(`and another polynomial in N. with coeffs. that are polynomials in x`):
 
fi:
 
 
 
 
 
 
if nops([args])=1 and op(1,[args])=`findrec2` then
print(` findrec2(f,DEGREEn,DEGREEk,ORDER,n,k,N)`):
print(` finds empirically an ordi. linear recurrence in n`):
print(` with polynomial coeffs. The input is a double sequence`):
print(`f[n,k] given as a list of  lists.`):
print(`STARTING at f[1,1],i.e. f[0] is not considered`):
print(` where DEGREEn:=the maximal degree of the coefficients in n`):
print(` where DEGREEk:=the maximal degree of the coefficients in k`):
print(`and ORDER:=the order of the recurrence in n.`):
print(`  The output is the operator`):
print(` in n, k and N, where N is the forward unit shift: Nf(n):=f(n+1).`):
print(`For example findrec2([[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]],0,0,1,n,k,N)`):
print(` should yield N-1`):
fi:
 
if nops([args])=1 and op(1,[args])=`Findrec` then
print(`Findrec(f,n,N,C) finds empirically an ordi. linear recurrence`):
print(` with polynomial coeffs. with ORDER+DEGREE<=C. The input is a sequence f given as a list`):
print(`STARTING at f[1],i.e. f[0] is not considered`):
print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`):
print(`For example Findrec([1,1,2,3,5,8,13,21,34],n,N) should yield`):
print(`N^2-N-1 , and findrec([1,2,5,14,42,132,429],m,M) should yield`):
print(`(m+1)*M-(4*m-2). If there is not enough data, you will get 0`):
print(`Of course, if a sequence is not P-recursive there would never`):
print(`enough data`):
fi:

if nops([args])=1 and op(1,[args])=`qfindrec` then
print(`qfindrec(f,DEGREE,ORDER,n,N) tries to finds empirically an ordi. linear recurrence`):
print(` with polynomial coeffs. in (q,q^n) a recurrence of degree DEGREE and order ORDER.`):
print(` The input is a sequence f given as a list`):
print(`of rational functions in q`):
print(`STARTING at f[1],i.e. f[0] is not considered`):
print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`):
print(`For example:`):
print(` qfindrec([seq(q^i,i=1..20)],0,1,n,N); `):
print(` should yield N-q  .`):
print(`If there is no recurrence of the specified order and degree, it returns FAIL`):
fi:
 
if nops([args])=1 and op(1,[args])=`decomp` then
print(`decomp(pol,x): given a polynomial in x, pol, decides whether it is`):
print(`the algebraic product of a quadratic polynomial pol1, and another pol`):
print(`pol2. If yes, it returns pol1 and pol2, if not, it returns 0`):
fi:
 
if nops([args])=1 and op(1,[args])=`pol_in_terms_of_p` then
print(`pol_in_terms_of_p(resh,x) :given the list of power sums [p_1,..,p_deg]`):
print(`finds the monic polynomial, in x, whose roots are the power sums of`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`newts` then
print(`newts(pol,x): gives a list of the power sums of the roots pol(x)`):
print(`from p_1 to p_(degree of pol(x)), using newtson's equations`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`newt` then
print(`newt(pol,x,R): gives a list of the power sums of the roots pol(x)`):
print(`from p_1 to p_R, using newtson's equations`):
fi:
 
if nops([args])=1 and op(1,[args])=`algprod` then
print(`algprod(P,Q,x): inputs polynomials P and Q in x , and outputs`):
print(`a poly R, in x, of degree deg(P)deg(Q) whose roots are all the`):
print(`possible products of the roots of P and Q`):
fi:
 
if nops([args])=1 and op(1,[args])=`logalg` then
 print(`logalg(F,P,x,ORDER): Using Comtet's method, given a polynomial`): 
 print(`F(P,x), finds the linear diff. eq. of order ORDER satisfied by the`): 
 print(`f.p.s. log(f(x)), where f(x) satisfies the alg.eq. F(f(x),x)=0 `):
 print(`E.g. try logalg(P-1-x*P^2,P,x,3) to get the differential `):
 print(`operator annihilating the log of the g.f for the Catalan numbers`):
 print(`ORDER=degree(F,P)+1 should always work, but you may be able to`):
 print(`to get away with lower order`):
fi:
 
if nops([args])=1 and op(1,[args])=`difftorec` then
 print(`difftorec(ope1,D,x,n,N): Converts the differential eq. (linear)`):
 print( `with poly coeff. satisfied by the g.f f(x) given by the opertor`):
 print(`ope (D,x),to the  recurrence, in n and N  satisfied `):
 print(`by the coeffs. of f(x). For example`):
 print(`difftorec(x*D-1,D,x,n,N) will yield n-1`):
fi:
 
if nops([args])=1 and op(1,[args])=`algtodiff` then
 print(`algtodiff(F,P,x,ORDER): Using Comtet's method, given a polynomial`): 
 print(`F(P,x), finds the linear diff. eq. of order ORDER satisfied by the`): 
 print(`formal power series f(x) that satisfies the alg.eq. F(f(x),x)=0 `):
 print(`Given an algebraic equation F(P(x),x) and desired`):
 print(`order, finds, if one exists, a linear diff. eq. of order ORDER`):
 print(`For example try algtodiff(P-1-x*P^2,P,x,2) to get the differential `):
 print(`operator annihilating the g.f for the Catalan numbers`):
 print(`if you then do difftorec(",x,D,n,N), you would get a new proof`):
 print(`that the Catalan numbers, defined by their combinatorial meaning`):
 print(`indeed equals (2n)!/(n!(n+1)!)`):
 print(`ORDER=degree(F,P) should always work, but you may be able to`):
 print(`to get away with lower order`):
 print(`Warning: it does not always gives the MINIMAL operator`):
fi:
 
if nops([args])=1 and op(1,[args])=`genseq` then
print(`genseq(eq,P,x,K) Generates the coefficients from x^0`):
print(`through x^K  of the formal`):
print(` powers series that is a solution of the eq P=eq(x,P)`):
print(`For example genseq(1/(1-x-x^2)-1,P,x,4) should yield`):
print(`[1,1,2,3,5] and genseq(1+x*P^2,P,x,5) should yield [1,1,2,5,14,42]`):
fi:
 
if nops([args])=1 and op(1,[args])=`genseqnew` then
print(`genseqnew(eq,P,x,K,PHI) outputs, the coeffs. from x^0 to x^K, as`):
print(`a list of polynomials in t (a global variable) of the f.p.s. `):
print(`PHI(x,t) that is a solution of the algebraic functional equation`):
 print(`PHI(x,t)=eq(PHI(x,t),P(x),x,t), where P(x)=PHI(x,1)`):
print(`Such a functional equation came up in Doron Zeilbeger's proof`):
print(`of Julian West's conj. about 2-Stack-sortable permutations.`):
print(`See Discrete Math 102(1992) 85-93. For example, try `):
print(` genseqnew(PHI-1/(1-x*t)-x*t*(P-t*PHI)*(P-PHI)/(1-t)^2,P,x,5,PHI)`):
print(`to get the first 6 terms in the f.p.s. of interest in the`): 
print(`above paper`):
print(`by taking K=40 and plugging in t=1, you can get a semi-rigorous`):
print(`proof of Julian West's conj. Of course the above paper has a `):
print(`proof`):
fi:
 
if nops([args])=1 and op(1,[args])=`empir` then
print(`empir(gu,degx,degP,x,P) empirically finds an algebraic equation`):
print(` F(P(x),x)=0 of degree degP in P(x) and degx in x for `):
print(`formal power series P(x):=sum_i gu[i]*x^i, where gu is a list`):
print(`for example empir([seq(1,i=1..20)],1,1,x,P) should yield`):
print(` P-xP-1=0`):
 
fi:
 
 
if nops([args])=1 and op(1,[args])=`Empir` then
print(`Empir(gu,x,P) empirically finds an algebraic equation`):
print(` F(P(x),x)=0 for the  `):
print(`formal power series P(x):=sum_{i=0} gu[i]*x^i, where gu is a list`):
print(`for example Empir([seq(1,i=1..20)],x,P) should yield`):
print(` P-xP-1=0`):
print(`If there is not enough data, it reurns 0. You should then`):
print(`try it again with a longer sequence`):
 
fi:
 
 
if nops([args])=1 and op(1,[args])=`mfindrec` then
print(`mfindrec(DEG,ORDER,f,p): like findrec(DEG,ORDER,f), but mod p`):
fi:
 
end:
 
pashet:=proc(p,N)
local i,gu1,gu,p1,ra:
p1:=normal(p):
gu1:=denom(p1):
ra:=degree(gu1,N):
p1:=subs(n=n+ra,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):
end:
 
 
 
difftorec:=proc(ope1,D,x,n,N)
local i,j,gu,ord,ope,r,mu:
 
ope:=expand(ope1):
gu:=0:
ord:=degree(ope,D):
 
for i from 0 to  ord do
 mu:=coeff(ope,D,i):
 
for j from 0 to degree(mu,x) do
  gu:=gu+coeff(mu,x,j)*product(n+r-j,r=1..i)*N^(i-j):
od:
od:
 
 
pashet(gu,N):
end:
 
 
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:=proc(F,P,x,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]*D:
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]*D^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):
ope:=expand(ope):
pip:=coeff(ope,D,degree(ope,D)):
 
pip:=coeff(pip,x,degree(pip,x)):
normal(ope/pip):
end:
 
 
 
 
genseq:=proc(eq,P,x,K)
local i,gu,P1:
gu:=[coeff(coeff(eq,P,0),x,0)]:
P1:=coeff(eq,P,0):
 
for i from 1 to K do
 P1:=subs(P=P1,eq):
 P1:=taylor(P1,x=0,i+1):
 gu:=[op(gu),coeff(P1,x,i)]:
od:
 
 
gu:
 
end:
 
   
 
 
#genseqnew:=proc(eq,P,x,K,PHI)
#PHI(x,t) that is a solution of the algebraic functional equation
#PHI(x,t)=eq(PHI(x,t),P(x),x,t), where P(x)=PHI(x,1)
 
genseqnew:=proc(eq,P,x,K,PHI)
local i,gu,P1,eq1,eq2,PHI1:
eq1:=normal(eq):
gu:=[1]:
PHI1:=1:
P1:=1:
for i from 1 to K do
 eq2:=subs(PHI=PHI1,eq1):
 eq2:=normal(eq2):
 eq2:=subs(P=P1,eq2):
 eq2:=normal(eq2):
PHI1:=taylor(eq2,x=0,i+1):
 PHI1:=normal(PHI1):
 P1:=subs(t=1,PHI1):
 
 gu:=[op(gu),coeff(PHI1,x,i)]:
od:
 
gu:
 
end:
 
 
 
#logalg(F,P,x,ORDER): Using 
#an adaptation of Comtet's method, given a polynomial
#F(P,x), finds the linear diff. eq. of order ORDER satisfied by the
#formal power series log (f(x)), where f(x) satisfies the alg.eq. F(f(x),x)=0 
 
 
 
logalg:=proc(F,P,x,ORDER)
local i,gu,a,lu,pa,deg,eq1,lu1,lu1a,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:=0:
lu:=subs(P=P(x),F):
lu1:=diff(lu,x):
lu1:=subs(diff(P(x),x)=y,lu1):
lu1:=solve(lu1=0,y):
lu1a:=lu1/P(x):
lu1a:=simp(pa,deg,lu1a,P,x):
 lu1a:=normal(lu1a):
 
 
gu:=gu+a[1]*lu1a:
gu:=normal(gu):
 gu:=simp(pa,deg,gu,P,x):
 
var:={a[1]}:
ope:=a[1]*D:
lu:=lu1a:
 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]*D^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):
ope:=expand(ope):
#pip:=coeff(ope,D,degree(ope,D)):
 
#pip:=coeff(pip,x,degree(pip,x)):
#normal(ope/pip):
end:
 
 
#algprod(P,Q,x): inputs polynomials P and Q in x , and outputs
#a poly R, in x, of degree deg(P)deg(Q) whose roots are all the
#possible products of the roots of P and Q
 
algprod:=proc(P,Q,x)
local degp,degq,deg,pol,var,eq,a,i,mu,lu,reshp,reshq,xp,xq,pol1,i1,i2:
 
degp:=degree(P,x):
 
degq:=degree(Q,x):
 
deg:=degp*degq:
 
reshp:=[]:
 
for i from 0 to degp-1 do
 reshp:=[op(reshp),x^i]:
od:
 
mu:=P/coeff(P,x,degp):
mu:=x^degp-mu:
 
reshp:=[op(reshp),mu]:
 
lu:=mu:
 
for i from 1 to deg-degp do
lu:=expand(x*lu):
lu:=subs(x^(degp)=mu,lu):
lu:=expand(lu):
reshp:=[op(reshp),lu]:
od:
 
 
reshq:=[]:
 
for i from 0 to degq-1 do
 reshq:=[op(reshq),x^i]:
od:
 
mu:=Q/coeff(Q,x,degq):
mu:=x^degq-mu:
 
reshq:=[op(reshq),mu]:
 
lu:=mu:
 
for i from 1 to deg-degq do
lu:=expand(x*lu):
lu:=subs(x^(degq)=mu,lu):
lu:=expand(lu):
reshq:=[op(reshq),lu]:
od:
 
reshp:=subs(x=xp,reshp):
reshq:=subs(x=xq,reshq):
pol:=0:
pol1:=0:
eq:={}:
var:={}:
 
for i from 0 to deg do
pol:=pol+a[i]*x^i:
var:=var union {a[i]}:
pol1:=expand(pol1+a[i]*op(i+1,reshp)*op(i+1,reshq)):
od:
 
 
for i1 from 0 to degp-1 do
for i2 from 0 to degq-1 do
 eq:=eq union {coeff(coeff(pol1,xp,i1),xq,i2)}:
od:
od:
 
var:=solve(eq,var):
pol:=subs(var,pol):
normal(pol/coeff(pol,x,degree(pol,x))):
end:
 
newts:=proc(pol,x)
local mu,deg,resh,pol1,r,i:
 
 
deg:=degree(pol,x):
 
pol1:=pol/coeff(pol,x,deg):
pol1:=expand(pol1):
resh:=[-coeff(pol1,x,deg-1)]:
 
 
for r from 2 to deg do
  mu:=0:
 
 for i from 1 to r-1 do
    mu:=mu-coeff(pol1,x,deg-i)*op(r-i,resh):
 od:
   
mu:=mu-r*coeff(pol1,x,deg-r):
mu:=expand(mu):
resh:=[op(resh),mu]:
 
od:
resh:
end:
 
 
 
newt:=proc(pol,x,R)
local mu,deg,resh,pol1,r,i:
 
 
if R<=degree(pol,x) then
  RETURN([op(1..R,newts(pol,x))]):
fi:
 
 
deg:=degree(pol,x):
 
pol1:=pol/coeff(pol,x,deg):
pol1:=expand(pol1):
resh:=newts(pol,x):
 
for r from deg+1 to R do
  mu:=0:
 
 for i from 1 to deg do
    mu:=mu-coeff(pol1,x,deg-i)*op(r-i,resh):
 od:
mu:=expand(mu):
resh:=[op(resh),mu]:
 
od:
resh:
end:
 
 
 
 
#pol_in_terms_of_p(resh) :given the list of power sums [p_1,..,p_deg]
#finds the monic polynomial whose roots are the power sums of
 
 
pol_in_terms_of_p:=proc(resh,x)
local pol,mu,deg,r,i:
 
deg:=nops(resh):
 
pol:=x^deg:
 
for r from 1 to deg do
  mu:=-op(r,resh):
 
 for i from 1 to r-1 do
    mu:=mu-coeff(pol,x,deg-i)*op(r-i,resh):
 od:
   
mu:=mu/r:
mu:=expand(mu):
pol:=pol+mu*x^(deg-r):
 
od:
pol:
end:
 
 
 
decomp:=proc(pol,x)
local i,deg,resh1,resh2,resh,a,b,pol1,pol2,eq,var,resh2a,eq1:
 
deg:=degree(pol,x):
 
deg:=deg/2:
 
if not type(deg, integer) then
  RETURN(0):
fi:
 
 
pol1:=x^2+a*x+b:
 
var:={a,b}:
 
resh1:=newt(pol1,x,2*deg):
resh:=newts(pol,x):
 
resh2:=[]:
 
for i from 1 to 2*deg do
 resh2:=[op(resh2),op(i,resh)/op(i,resh1)]:
od:
 
pol2:=pol_in_terms_of_p([op(1..deg,resh2)],x):
 
resh2a:=newt(pol2,x,2*deg):
 
eq:={}:
 
for i from 1 to deg do
 eq1:=normal(op(i,resh2)-op(i,resh2a)):
  if eq1<>0 then
    ERROR(`Something is wrong`):
  fi:
od:
 
 
for i from deg+1 to 2*deg do
 eq1:=normal(op(i,resh2)-op(i,resh2a)):
  eq:=eq union {eq1}:
od:
 
var:=solve(eq,var):
 
if var=NULL then
 RETURN(0):
fi:
 
pol1:=subs(var,pol1):
 
pol2:=subs(var,pol2):
 
pol1,pol2:
 
end:
 
 
decompx:=proc(pol,N,x)
local i1,i,mu,deg,resh1,resh2,resh,a,pol1,pol2,eq,var,resh2a,eq1:
 
deg:=degree(pol,N):
 
deg:=deg/2:
 
if not type(deg, integer) then
  RETURN(0):
fi:
 
 
pol1:=N^2+2*a*x*N-1:
 
var:={a}:
 
resh1:=newt(pol1,N,2*deg):
resh:=newts(pol,N):
 
resh2:=[]:
 
for i from 1 to 2*deg do
 resh2:=[op(resh2), rem(op(i,resh),op(i,resh1),x)]:
od:
 
eq:={}:
 
for i from 1 to 3 do
 mu:=op(i,resh2):
  mu:=expand(numer(normal(mu))):
 for i1 from 0 to degree(mu,x) do
  eq:=eq union {numer(normal(coeff(mu,x,i1)))=0}:
 od:
od:
var:=solve(eq,var):
 
if var=NULL then
 RETURN(0):
fi:
 
var:
 
end:
 
 
 
 
findrec2:=proc(f,DEGREEn,DEGREEk,ORDER,n,k,N)
local ope,var,eq,a,i,i1,j,n0,k0,kv,var1,eq1,mu:
if (1+DEGREEn)*(1+DEGREEk)*(2+ORDER)+1+ORDER>nops(f)*nops(op(1,f)) then
ERROR(` Not enough date`):
fi:
 
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for i1 from 0 to DEGREEn do
 for j from 0 to DEGREEk do
  ope:=ope+a[i,i1,j]*n^i1*k^j*N^i:
  var:=var union {a[i,i1,j]}:
 od:
od:
od:    
eq:={}:
 
mu:=trunc(evalf(sqrt((1+DEGREEn)*(1+ORDER)*(1+DEGREEk)))):
 
#if mu+2+ORDER
for n0 from 1 to mu+2 do
for k0 from 1 to mu+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs({k=k0,n=n0},coeff(ope,N,i))*op(k0,op(n0+i,f)):
  od:
 
   eq:= eq union {eq1}:
od:
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 nops(kv)>1 then
  print(` either DEGREE or ORDER are too high`):
   print(`The output is not the minimal possible operator`):
fi:
 
for i from 1 to nops(kv) do
  ope:=subs(op(i,kv)=1,ope):
od:
 
ope:
 
end:
 
 
###############From Findrec 

with(linalg):
ka:=proc():
if args=0 then
print(`kaki`):

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);`):
else
print(`Hi`):
fi:
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)]:
od:

gu:

end:


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

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

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:

#KamaArokh(C): the length of the sequence needed to guess a recurrence of complexity C
KamaArokh:=proc(C) local DEGREE:
max(seq((2+DEGREE)*(1+C-DEGREE)+5,DEGREE=0..C)):
end:


 
###############End from Findrec
 
 
 


#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,flo,pip,var2,mu1,halev,vu,i,vu1:
if (1+degx)*(1+degP) > nops(gu)-3 then
 RETURN(`sequence too small`):
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,nops(gu)-1):
 
eq:={}:
 
for i1 from 0 to nops(gu)-2 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):
 


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

var2:={}:

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

if nops(var2)<>1 then
F:=coeff(F,var2[1],1):
fi:



flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=F/pip:
F:=numer(normal(F)):

vu:=add(gu[i+1]*x^i,i=0..nops(gu)-1):

halev:=add(factor(coeff(F,P,i))*P^i,i=0..degree(F,P)):

vu1:=taylor(subs(P=vu,halev),x=0,nops(gu)+1):

if {seq(coeff(vu1,x,i),i=1..nops(gu)-1)}<>{0} then
RETURN(0):
fi:

halev:

end:

Empir:=proc(gu,x,P)
local degx,degP,L,lu:
 
for L from 1 to (nops(gu)-3)/3 do
for degP from 1 to L do
for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do
 
lu:=empir(gu,degx,degP,x,P):
 
if lu<>0 then
RETURN(lu):
fi:
od:
od:
od:
0:
end:
 









#######End from SCHUTZENBERGER

#OneStep(F,x,z): given a formal power series F in x,z, performs one iteration in
#F->x+x*T(F), where T(f)=(subs(z=1,f)-z^2*f)/(1-z)
#Try:
#OneStep(x,x,z);
OneStep:=proc(F,x,z) local gu:
gu:=normal(F/(1-F)):
normal(x+x*(subs(z=1,gu)-z^2*gu)/(1-z)):

end:


#OneStepK(F,x,z,K): given a formal power series F in x,z, performs one iteration in
#F->x+x*T(F), where T(f)=(subs(z=1,f)-z^2*f)/(1-z) and truncates it to the K-th term in z
#Try:
#OneStepK(x,x,z,K);
OneStepK:=proc(F,x,z,K) local gu,i:
gu:=normal(F/(1-F)):
gu:=taylor(gu,x=0,K+1):
gu:=add(expand(coeff(gu,x,i))*x^i,i=0..K):
gu:=normal(x+x*(subs(z=1,gu)-z^2*gu)/(1-z)):
gu:=add(expand(coeff(gu,x,i))*x^i,i=0..K):
gu:
end:


#PreT01(K,z): The first K terms in the sequence of the generating function, according to the label of the root, of Pre-Beta01 trees
#Try:
#PreT01(10,z);
PreT01:=proc(K,z) local gu,i,x:
gu:=x:


for i from 1 to K do
 gu:=OneStepK(gu,x,z,K):
od:

[seq(coeff(gu,x,i),i=1..K)]:

end:

#T01(K,z): The first K terms in the sequence of the generating function, according to the label of the root, of Beta01 trees
#Try:
#T01(10,z);
T01:=proc(K,z) local gu,i,x:
gu:=x:


for i from 1 to K do
 gu:=OneStepK(gu,x,z,K):
od:

gu:=x/(1-gu):

gu:=taylor(gu,x=0,K+1):


[seq(expand(coeff(gu,x,i)),i=1..K)]:

end:


#FunEq01(F,x,z): The Functional Equation for F(x,z): the genrating function for Pre-01-trees where the weight is x^(#vertices)*z^(label of root), Try:
#FunEq01(F,x,z)
FunEq01:=proc(F,x,z) local gu:

gu:=(F(x,1)/(1-F(x,1))-z^2*F(x,z)/(1-F(x,z)))/(1-z):

F(x,z)=x+x*gu:

end:



#CheckFunEq01(F,x,z,K): Checks that F satisfies  the Functional Equation for F(x,z) for Pre-01-trees up to K terms, Try:
#gu:=PreT01(20,z):
#gu:=add(gu[i]*x^i,i=1..20):
#CheckFunEq01(gu,x,z,20);
CheckFunEq01:=proc(F,x,z,K) local gu,i:

gu:=normal(F/(1-F)):

gu:=normal(normal(x+x*(subs(z=1,gu)-z^2*gu)/(1-z))-F):

gu:=taylor(gu,x=0,K+1):

evalb({seq(expand(coeff(gu,x,i)),i=1..K)}={0}):


end:






#PreBetaEq01(P,x,z,K) : Guessing an algebraic equation satisfied by the generating function of Pre-Beta-01-trees, using K terms, try:
#PreBetaEq01(P,x,z,40);
PreBetaEq01:=proc(P,x,z,K)  local gu,gu2,gu3,eq2,eq3,degx,degP,eq:

gu:= PreT01(K,z):

gu2:=subs(z=2,gu):

eq2:=Empir(gu2,x,P):

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

degx:=degree(eq2,x):
degP:=degree(eq2,P):

gu3:=subs(z=3,gu):
eq3:=Empir(gu3,x,P):

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

if degree(eq3,x)<>degx or degree(eq3,P)<>degP then
print(` Inconsistend degrees `):
RETURN(FAIL):
fi:


eq:=empir(gu,degx,degP,x,P):

if eq=0 then
 RETURN(FAIL):
else
RETURN(eq):
fi:

end:

#BetaEq01(P,x,z,K) : Guessing an algebraic equation satisfied by the generating function of Pre-Beta-01-trees, using K terms, try:
#BetaEq01(P,x,z,40);
BetaEq01:=proc(P,x,z,K)  local gu,gu2,gu3,eq2,eq3,degx,degP,eq:

gu:= T01(K,z):

gu2:=subs(z=2,gu):

eq2:=Empir(gu2,x,P):

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

degx:=degree(eq2,x):
degP:=degree(eq2,P):

gu3:=subs(z=3,gu):
eq3:=Empir(gu3,x,P):

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

if degree(eq3,x)<>degx or degree(eq3,P)<>degP then
print(` Inconsistend degrees `):
RETURN(FAIL):
fi:


eq:=empir(gu,degx,degP,x,P):

if eq=0 then
 RETURN(FAIL):
else
RETURN(eq):
fi:

end:



 
 
 






#AsyAnal(L,N,n): given a list of numbers L, of length k, say, whose i-th entry is f(N+i-1) for i from 1 to nops(L), conjectures
#an asymptotic approximate expression for f(m) in the form a[0]+a[1]/n+...+ a[k]/n^(k-1). Try:
#AsyAnal([seq(1+6/i+11/i^2,i=98..100)],98,n);
AsyAnal:=proc(L,N,n) local k,eq,var,X,a,i:

k:=nops(L)-1:

X:=add(a[i]/n^i,i=0..k):

var:={seq(a[i],i=0..k)}:

eq:={seq(subs(n=N+i-1,X)-L[i],i=1..k+1)}:


var:=solve(eq,var):

subs(var,X):

end:



#OneStepBab(F,x,y,a,b): given a formal power series F in x,y, performs one iteration in
#the evolution of Beta(a,b)-Description trees according to the article:
#F->x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F)
#Following the article:
#Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer. 
#Description trees for some families of planar maps." Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.
#http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz
#Try:
#OneStepBab(x,x,z,0,1);
OneStepBab:=proc(F,x,y,a,b):
normal(x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F))):

end:

#OneStepBabK(F,x,y,a,b,K): given a formal power series F in x,y, performs one iteration, and truncates it to K terms in x, in
#the evolution of Beta(a,b)-Description trees 
#F->x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F)
#Following the article:
#Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer. 
#Description trees for some families of planar maps." Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.
#http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz
#Try:
#OneStepBabK(x,x,z,0,1,10,10);
OneStepBabK:=proc(F,x,y,a,b,K) local gu,i:
gu:=normal(x*(1+F)*(y^a+1/(1-y)*(y^a*subs(y=1,F)-y^(b+1)*F))):
gu:=taylor(gu,x=0,K+1):
gu:=add(expand(coeff(gu,x,i))*x^i,i=0..K):
end:




#Bab(a,b,y,K): The first K terms in the sequence of the generating function, according to the label of the root, of 
#Beta(a,b) Description trees try:
#Bab(0,1,z,10);
Bab:=proc(a,b,y,K) local gu,i,x:
option remember:
gu:=0:


for i from 1 to K do
 gu:=OneStepBabK(gu,x,y,a,b,K):
od:

[seq(coeff(gu,x,i),i=1..K)]:

end:



#BabFE(a,b,x,y,fxy,fx1): the functional equation for B(a,b) description trees, given
#as a function of x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1)),
#the functional equation itself is
#f(x,y)=F(x,y,f(x,y),f(x,1)); Try:
#BabFE(0,1,x,y,fxy,fx1);
BabFE:=proc(a,b,x,y,fxy,fx1)
x*(1+fxy)*(y^a+1/(1-y)*(y^a*fx1-y^(b+1)*fxy)):
end:

#AabFE(a,b,x,y,fxy,fx1): the functional equation for A(a,b) description trees, given
#as a function of x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1)),
#the functional equation itself is
#f(x,y)=F(x,y,f(x,y),f(x,1)); Try:
#AabFE(0,1,x,y,fxy,fx1);
AabFE:=proc(a,b,x,y,fxy,fx1)
x*y^b+1/(1-fxy)*x/(1-y)*(y^a*fx1-y^(b+1)*fxy):
end:


#jwFE(x,y,fxy,fx1): The Functional Equation for Two-Stack-Sortable Permutations
jwFE:=proc(x,y,fxy,fx1):
1/(1-x*y)+x*y*(fx1-y*fxy)*(fx1-fxy)/(1-y)^2:
end:

#GenjwFE(x,y,fxy,fx1,A,B): The Generalization of the Functional Equation for Two-Stack-Sortable Permutations
#with two extra parameters, A,B. When A=1, B=1 it is the original. Try:
#GenjwFE(x,y,fxy,fx1,2,3);
GenjwFE:=proc(x,y,fxy,fx1,A,B):
A/(1-x*y)+B*x*y*(fx1-y*fxy)*(fx1-fxy)/(1-y)^2:
end:


#FirstK(F,x,y,fxy,fx1,K): The list of coeficients of x^i, from i=1 to i=K 
#of the formal power series in x,y, that satisfies the functional equation
#f(x,y)=F(x,y,f(x,y),f(x,1)). Try:
#FirstK(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,10);
FirstK:=proc(F,x,y,fxy,fx1,K)
local gu,i:
option remember:
gu:=0:


for i from 1 to K do
 gu:=OneStepFeK(gu,F,x,y,fxy,fx1,K):
od:

[seq(coeff(gu,x,i),i=0..K)]:

end:

#OneStepFeK(g,F,x,y,fxy,fx1,K): given a formal power series g applies one 
#Try:
#OneStepFeK(0,BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,10);
OneStepFeK:=proc(g,F,x,y,fxy,fx1,K) local gu,i:

gu:=normal(subs({fxy=g,fx1=subs(y=1,g)},F)):


gu:=taylor(gu,x=0,K+1):
gu:=add(expand(coeff(gu,x,i))*x^i,i=0..K):

end:


#FEtoAE1(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))
#guesses, using the first K terms, an algebraic equation for f(x,1). Try:
#FEtoAE1(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);
FEtoAE1:=proc(F,x,y,fxy,fx1,K) local gu,ka:
option remember:
gu:=subs(y=1,FirstK(F,x,y,fxy,fx1,K+10)):

gu:=[op(1..K,gu)]:

if gu[nops(gu)]=0 then
 gu:=[op(1..nops(gu)-1,gu)]:
fi:

if gu[nops(gu)]=0 then
 gu:=[op(1..nops(gu)-1,gu)]:
fi:

if gu[nops(gu)]=0 then
 gu:=[op(1..nops(gu)-1,gu)]:
fi:



ka:=Empir(gu,x,fx1):


if ka=0 then
 RETURN(FAIL):
else
 RETURN(ka):
fi:

end:

#FEtoRec1(F,x,y,fxy,fx1,n,N,C): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))
#guesses a recurrence operator in terms of n and the shift operator N of comlexity <=C annihilating the
#coefficents of f(x,1). Try:
#FEtoRec1(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,n,N,6);
FEtoRec1:=proc(F,x,y,fxy,fx1,n,N,C) local gu,K,D1:

K:=max(seq((1+D1)*(1+C-D1)+2+C-D1,D1=0..C)):
gu:=subs(y=1,FirstK(F,x,y,fxy,fx1,K)):


Findrec([op(2..nops(gu),gu)], n,N,C):

end:

#FEtoAEslow(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))
#guesses, using the first K terms, an algebraic equation for f(x,y). Try:
#FEtoAEslow(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);
FEtoAEslow:=proc(F,x,y,fxy,fx1,K) local gu, gu2,gu3,dx,dP,P:

gu:=FirstK(F,x,y,fxy,fx1,K):

gu2:=Empir(subs(y=2,gu),x,P):

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

dx:=degree(gu2,x):
dP:=degree(gu2,P):

gu3:=Empir(subs(y=3,gu),x,P):

if (degree(gu3,x)<>dx or degree(gu3,P)<>dP) then
RETURN(FAIL):
fi:


empir(gu,dx,dP,x,fxy):

end:


#FEtoAE(F,x,y,fxy,fx1,K): Using Ira Gessel's idea.
#Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))
#guesses, using the first K terms, an algebraic equation for f(x,y). Try:
#FEtoAE(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);
FEtoAE:=proc(F,x,y,fxy,fx1,K) local A,B,gu,mu,i,lu,a0:
option remember:
A:=numer(normal(fxy-F)):
B:=FEtoAE1(F,x,y,fxy,fx1,K):


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

gu:=expand(resultant(A,B,fx1)):


mu:=gu:
while subs(x=0,mu)=0 do
 mu:=normal(mu/x):
od:
gu:=mu:


mu:=FirstK(F,x,y,fxy,fx1,20):


a0:=mu[1]:
mu:=add(mu[i+1]*x^i,i=0..19):

mu:=expand(taylor(subs(fxy=mu,gu),x=0,20)):




if {seq(expand(coeff(mu,x,i)),i=0..18)}<>{0} then
print(`Something terrible happend`):
RETURN(FAIL):
else
RETURN(gu):
fi:



end:




#FEtoAE1wp(F,x,y,fxy,fx1,K): Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))
#guesses  using the first K terms, an algebraic equation for f(x,1). It then proceeds to prove it!
#If it fails to prove it it says so.
#Try:
#FEtoAE1wp(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30);
FEtoAE1wp:=proc(F,x,y,fxy,fx1,K) local gu,lu:
option remember:
gu:=subs(y=1,FirstK(F,x,y,fxy,fx1,K+10)):

gu:=[op(1..K,gu)]:

if gu[nops(gu)]=0 then
 gu:=[op(1..nops(gu)-1,gu)]:
fi:

if gu[nops(gu)]=0 then
 gu:=[op(1..nops(gu)-1,gu)]:
fi:

if gu[nops(gu)]=0 then
 gu:=[op(1..nops(gu)-1,gu)]:
fi:


gu:=Empir(gu,x,fx1):

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

lu:=FEtoAE(F,x,y,fxy,fx1,K):

lu:=subs({fxy=fx1,y=1},lu):

lu:=denom(normal(lu/gu)):

if degree(lu,fx1)>0 then
print(gu):
print(`did not work out`):
RETURN(FAIL):
else
RETURN(gu):
fi:


end:


#PrintFE(F,x,y,fxy,fx1,f) : prints the functional equation nicely
PrintFE:=proc(F,x,y,fxy,fx1,f) :
print(f(x,y)=subs({fxy=f(x,y),fx1=f(x,1)},F)):
end:

#FEtoAE1wpV(F,x,y,fxy,fx1,K,f): Verbose version of FEtoAE(F,x,y,fxy,fx1,K). Using Ira Gessel's idea.
#Given a functional equation, F, an expression in x,y,fxy (denoting f(x,y)) and fx1 (denoting f(x,1))
#guesses, using the first K terms, an algebraic equation for f(x,y). Try:
#FEtoAE1wpV(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30,f);
FEtoAE1wpV:=proc(F,x,y,fxy,fx1,K,f) local A,B,gu,mu,i,lu,t0:
t0:=time():

A:=numer(normal(fxy-F)):

B:=FEtoAE1(F,x,y,fxy,fx1,K):

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

gu:=factor(resultant(A,B,fx1)):

mu:=FirstK(F,x,y,fxy,fx1,8):

mu:=add(mu[i+1]*x^i,i=0..7):

mu:=expand(subs(fxy=mu,gu)):

if ldegree(mu,x)<7 then
print(`Something terrible happend`):
RETURN(FAIL):
fi:

print(` Solution of a Certain Functional Equation`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let`, f(x,y), ` be the formal power series, in`, x,y, `satisfying the FUNCTIONAL equation`):
PrintFE(F,x,y,fxy,fx1,f):
print(``):
print(`then P=f(x,1) satisfies the ALGEBRAIC equation`):
print(subs(fx1=f(x,1),B)=0):
print(``):
print(`Proof: `):
print(`Moving everything in the functional equation to the left, clearing denominators and taking the numerator entails the followsing algebraic relation`):
print(`between `, f(x,y), f(x,1),x, ` and`,  y ):
print(subs({fxy=f(x,y),fx1=f(x,1)},A)=0):

gu:=factor(resultant(A,B,fx1)):

print(`Eliminating f(x,1) between these two equations, by taking the resultant, gives the following ALGEBRAIC relation between `, f(x,y),x,y ):

print(subs(fxy=f(x,y),gu)=0):

mu:=FirstK(F,x,y,fxy,fx1,8):

mu:=add(mu[i+1]*x^i,i=0..7):

mu:=expand(subs(fxy=mu,gu)):

if ldegree(mu,x)<7 then
print(`Something terrible happend`):
RETURN(FAIL):
fi:

print(`To prove that the above algebraic equation for f(x,1) , namely`):

print(subs(fx1=f(x,1),B)=0):

print(`is the right one, let's substitute`, y=1 ,`in the above relation for`, f(x,y), ` getting `):
gu:=subs({fxy=fx1,y=1},gu):
print(factor(subs(fx1=f(x,1),gu))=0):


lu:=FEtoAE(F,x,y,fxy,fx1,K):

lu:=subs({fxy=fx1,y=1},lu):

lu:=normal(lu/B):
print(`The quotient of the right sides is`, lu):

lu:=denom(lu):
if degree(lu,fx1)>0 then
print(gu):
print(`did not work out`):
RETURN(FAIL):
else
print(`the latter is indeed a multiple of the former, QED!`):
print(`This ends this article, that took`, time()-t0, `seconds. to generate. Have a great day!`):
fi:

end:







#########From Findrec
 








#AEtoSeq(F,P,x,a0,K): the first K coefficients of the solution of F(P,x)=0 with constant term a0, try:
#AEtoSeq(C-1-x*C^2,C,x,1,10);
AEtoSeq:=proc(F,P,x,a0,K) local gu,i,a,P1,a1,ka:

gu:=[]:
P1:=a0:

for i from 1 to K do
 P1:=P1+a*x^i:
 ka:=numer(normal(expand(coeff(taylor(subs(P=P1,F),x=0,i+1),x,i)))):


 a1:=solve({ka},{a}):
  if a1=NULL then

    RETURN(FAIL):
 else
   a1:=subs(a1,a):
  fi:
 gu:=[op(gu),a1]:
 P1:=subs(a=a1,P1):
od:

gu:

end:

 



#Info(F,x,y,fxy,fx1,K,n,N,MaxC,L0,L1): 
#Inputs:
#(i) A polynomial F in the variables x,y,fxy,fx1, representing a functional equation F(f(x,y),f(x,1),x,y)=0
#(ii) variables x,y,fx1,fx1, for expressing the above-mentioned Functional Equation
#(iii) a positive integer K indicating how many terms to use for guessing the algebraic equation
#(iv) symbols n and N for guessing a linear recurrence operator ope(n,N) annihilating the sequence
#(v): MaxC: a positive integer MaxC indicating the maximum complexity of the recurrence
#(vi) positive integers L0 and L1, 
#outputs (in the best possible case) a list [eq,ope,ListOfFirstL0 terms, L1-th term]
#Try:
#Info(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30,n,N,10,30,1000);

Info:=proc(F,x,y,fxy,fx1,K,n,N,C,L0,L1) local eq,a0,lu,K1,D1,gu,ope,ku:
option remember:

eq:=FEtoAE1(F,x,y,fxy,fx1,K):



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

a0:=subs(y=1,FirstK(F,x,y,fxy,fx1,1)[1]):

lu:=AEtoSeq(eq,fx1,x,a0,L0):

if lu=FAIL then
lu:=subs(y=1,FirstK(F,x,y,fxy,fx1,L0)):
RETURN([eq,lu]):
fi:

K1:=max(seq((1+D1)*(1+C-D1)+7+C-D1,D1=0..C)):

gu:=AEtoSeq(eq,fx1,x,a0,K1):


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


if ope=FAIL then
RETURN([eq,lu]):
fi:



if degree(ope,N)>0 then
if {seq(add(subs(n=n1,coeff(ope,N,i1))*gu[n1+i1],i1=0..degree(ope,N)),n1=5..nops(gu)-degree(ope,N))}<>{0} then
RETURN([eq,lu]):
fi:
fi:

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

[eq,lu,ope,ku]:

end:


#InfoV(F,x,y,fxy,fx1,K,n,N,MaxC,L0,L1): 
#verbose version of Info(F,x,y,fxy,fx1,K,n,N,MaxC,L0,L1)(q.v.)
#Try:
#InfoV(BabFE(0,1,x,y,fxy,fx1),x,y,fxy,fx1,30,n,N,10,30,1000);

InfoV:=proc(F,x,y,fxy,fx1,K,n,N,C,L0,L1) local gu,f,P,eq1,A,i,ope:

gu:=Info(F,x,y,fxy,fx1,K,n,N,C,L0,L1):

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

print(``):
print(`Theorem: Let f(x,y) be (unique!) formal power series, in the variables, x,y, satisfying the FUNCTIONAL Equation`):
PrintFE(F,x,y,fxy,fx1,f) :
eq1:=subs(fx1=P,gu[1]):
print(`Then P(x)=f(x,1) satisfies the following algebraic equation `):
print(``):
print(eq1=0):
print(``):
print(`and in Maple input notation`):
lprint(eq1=0):
print(``):
print(`The first`, L0, `terms, for the sake of the OEIS are `):
print(``):
lprint(gu[2]):
print(``):

if nops(gu)=2 then
RETURN(gu):
fi:

ope:=gu[3]:
print(`Furthermore, the sequence of coefficients, let's call them `,  A[n], `satisfy: `):
print(``):
print(add(coeff(ope,N,i)*A[n+i],i=0..degree(ope,N))=0):
print(``):
print(`and in Maple input format`):
print(``):
lprint(add(coeff(ope,N,i)*A[n+i],i=0..degree(ope,N))=0):
print(``):
print(`Just for fun`, A[L1], `equals `):
print(``):
lprint(gu[4]):
print(``):
RETURN(gu):
end:


#SeferAlpha(K,C,L0,L1,g): inputs positive integers K,C,L0,L1 and g and outputs an article about
#the enumeration of  Alpha(a,b) Description trees for 0<=a,b<=g. For each it also
#outputs the linear recurrence equation with polynomial coefficients satisfied by the enumerating sequence
#if its ORDER+DEGREE<=C.
#It also outputs the set of all  successful
#pairs [a,b] for which K data points suffice. For example, Try:
#SeferAlpha(20,10,30,100,2);

SeferAlpha:=proc(K,C,L0,L1,g) local gu,mu,t0,fxy,fx1,x,y,f,a,b,mu1,n,N:

t0:=time():

mu:=0:

print(`Enumeration of Alpha(a,b) Description Trees for  a and b up to`, g):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In this article, we will try to enumerate Alpha(a,b) description trees for a,b up to `, g):
print(`by discovering (easily provable, once discovered!) algebraic equations for the generating functions`):
print(`of the enumerating sequence, and will also present, in each case, the implied linear recurrence equation with polynomial coefficients`):
print(`satisfied by the enumerating sequence in each case, provided its order + degree are <=`, C):
print(`We will do this by using as TRAINING DATA (so to speak!) the first `, K, `terms, and admit failure otherwise`):
print(`For the sake of the OEIS, we will give the first`, L0, `terms in each case, and just for fun the`, L1, `-th term. `):
print(`Finally we will output the set of all successful pairs [a,b]. `):
print(`According to the article `):
print(`Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer. "Description trees for some families of planar maps" `):
print(`Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.`):
print(`Available from  http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz`):
print(`[viewed Dec. 24, 2014]`):
print(`The generating function, f(x,y),in the two variables (y being the so-called catalytic variable taking care of the label of the root)`):
print(`satisfies the FUNCTIONAL EQUATION`):
PrintFE(AabFE(a,b,x,y,fxy,fx1),x,y,fxy,fx1,f); 
print(`we would be interested in the straight enumeration, whose generating function is f(x,1) `):

mu:={}:
for a from 0 to g do
for b from 0 to g do


gu:=Info( AabFE(a,b,x,y,fxy,fx1) ,x,y,fxy,fx1,K,n,N,C,L0,L1):

if gu<>FAIL then
print(``):
print(`--------------------------------------------------------`):
print(``):
mu:=mu union {[a,b]}:

print(`For the Enumeration of `, Alpha(a,b), `trees we have the following `):
print(``):
InfoV( AabFE(a,b,x,y,fxy,fx1),x,y,fxy,fx1,K,n,N,C,L0,L1):
print(``):
fi:
od:
od:

mu1:={seq(seq([a,b],a=0..g),b=0..g)} minus mu:
print(``):
print(`----------------------------------------------------`):
print(``):
print(`This concludes this article for enuermationg Alpha(a,b) trees with a and b <=`,g):
print(`We were successful for the following pairs`, mu):

if mu1<>{} then
print(`but failed for the following pairs`, mu1):
else
print(`In other words, they all worked!`):
fi:
print(``):
print(`This colcludes this article, that took`, time()-t0, `seconds to generate. `):


end:


#SeferBeta(K,C,L0,L1,g): inputs positive integers K,C,L0,L1 and g and outputs an article about
#the enumeration of  Beta(a,b) Description trees for 0<=a,b<=g. For each it also
#outputs the linear recurrence equation with polynomial coefficients satisfied by the enumerating sequence
#if its ORDER+DEGREE<=C.
#It also outputs the set of all  successful
#pairs [a,b] for which K data points suffice. For example, Try:
#SeferBeta(20,10,30,100,2);

SeferBeta:=proc(K,C,L0,L1,g) local gu,mu,t0,fxy,fx1,x,y,f,a,b,mu1,n,N:

t0:=time():

mu:=0:

print(`Enumeration of Beta(a,b) Description Trees for  a and b up to`, g):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In this article, we will try to enumerate Beta(a,b) description trees for a,b up to `, g):
print(`by discovering (easily provable, once discovered!) algebraic equations for the generating functions`):
print(`of the enumerating sequence, and will also present, in each case, the implied linear recurrence equation with polynomial coefficients`):
print(`satisfied by the enumerating sequence in each case, provided its order + degree are <=`, C):
print(`We will do this by using as TRAINING DATA (so to speak!) the first `, K, `terms, and admit failure otherwise`):
print(`For the sake of the OEIS, we will give the first`, L0, `terms in each case, and just for fun the`, L1, `-th term. `):
print(`Finally we will output the set of all successful pairs [a,b]. `):
print(`According to the article `):
print(`Cori, Robert, Benjamin Jacquard, and Gilles Schaeffer. "Description trees for some families of planar maps" `):
print(`Proceedings of the 9th Conference on Formal Power Series and Algebraic Combinatorics. 1997.`):
print(`Available from  http://www-igm.univ-mlv.fr/~fpsac/FPSAC97/ARTICLES/Schaeffer.ps.gz`):
print(`[viewed Dec. 24, 2014]`):
print(`The generating function, f(x,y),in the two variables (y being the so-called catalytic variable taking care of the label of the root)`):
print(`satisfies the FUNCTIONAL EQUATION`):
PrintFE(BabFE(a,b,x,y,fxy,fx1),x,y,fxy,fx1,f); 
print(`we would be interested in the straight enumeration, whose generating function is f(x,1) `):

mu:={}:

for a from 0 to g do
for b from 0 to g do


gu:=Info( BabFE(a,b,x,y,fxy,fx1) ,x,y,fxy,fx1,K,n,N,C,L0,L1):

if gu<>FAIL then
print(``):
print(`--------------------------------------------------------`):
print(``):
mu:=mu union {[a,b]}:

print(`For the Enumeration of`,  BETA(a,b), `trees we have the following `):
print(``):
InfoV( BabFE(a,b,x,y,fxy,fx1),x,y,fxy,fx1,K,n,N,C,L0,L1):
print(``):
fi:
od:
od:

mu1:={seq(seq([a,b],a=0..g),b=0..g)} minus mu:
print(``):
print(`----------------------------------------------------`):
print(``):
print(`This concludes this article for enuermationg Beta(a,b) trees with a and b <=`,g):
print(`We were successful for the following pairs`, mu):

if mu1<>{} then
print(`but failed for the following pairs`, mu1):
else
print(`In other words, they all worked!`):
fi:
print(``):
print(`This colcludes this article, that took`, time()-t0, `seconds to generate. `):


end:



#SeferJW(K,C,L0,L1): inputs positive integers K,C,L0,L1 and g and outputs an article about
#the enumeration of  2-stack-sortable permutations
#Try
#SeferJW(20,10,30,100);

SeferJW:=proc(K,C,L0,L1) local gu,t0,fxy,fx1,x,y,f,n,N:

t0:=time():


print(`Enumeration of 2-Stack-Sortable Permutations `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In this article, we will  enumerate  2-Stack-Sortable Permutations`):
print(`by discovering (easily provable, once discovered!) an algebraic equation for the generating function`):
print(`of the enumerating sequence, and will also present the implied linear recurrence equation with polynomial coefficients`):
print(`satisfied by the enumerating sequence`):
print(`We will do this by using as TRAINING DATA (so to speak!) the first `, K, `terms, and admit failure otherwise`):
print(`According to the article `):
print(`Doron Zeilberger, A Proof of Julian West's Conjecture ... `):
print(`Discrete Math 102(1992), 85-93`):
print(`Available from   http://www.math.rutgers.edu/~zeilberg/mamarimY/julian.pdf `):
print(`The generating function, f(x,y),in the two variables (y being the so-called catalytic variable)`):
print(`satisfies the FUNCTIONAL EQUATION`):
PrintFE(jwFE(x,y,fxy,fx1),x,y,fxy,fx1,f); 
print(`we would be interested in the straight enumeration, whose generating function is f(x,1) `):


gu:=Info( jwFE(x,y,fxy,fx1) ,x,y,fxy,fx1,K,n,N,C,L0,L1):

if gu<>FAIL then
print(``):
print(`--------------------------------------------------------`):
print(``):

print(`For the Enumeration of 2-stack-sortable permutations, we have `):
print(``):
InfoV( jwFE(x,y,fxy,fx1),x,y,fxy,fx1,K,n,N,C,L0,L1):
print(``):
fi:


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

print(`This colcludes this article, that took`, time()-t0, `seconds to generate. `):


end:




#GenBabFE(a,b,A,B,C,x,y,fxy,fx1), a generalization of BabFE(a,b,x,y,fxy,fx1) with parameters A,B,C
#GenBabFE(a,b,1,1,1,x,y,fxy,fx1) is the same as BabFE(a,b,xy,fxy,fx1): Try:
#GenBabFE(0,1,1,1,1,x,y,fxy,fx1);
GenBabFE:=proc(a,b,A,B,C,x,y,fxy,fx1):
x*(A+fxy)*(B*y^a+C/(1-y)*(y^a*fx1-y^(b+1)*fxy)):
end:


#GenAabFE(a,b,A,B,x,y,fxy,fx1), a generalization of AabFE(a,b,x,y,fxy,fx1) with parameters A,B,C
#GenAabFE(a,b,1,1,x,y,fxy,fx1) is the same as AabFE(a,b,xy,fxy,fx1): Try:
#GenAabFE(0,1,1,1,x,y,fxy,fx1);
GenAabFE:=proc(a,b,A,B,x,y,fxy,fx1)
A*x*y^b+B/(1-fxy)*x/(1-y)*(y^a*fx1-y^(b+1)*fxy):
end:
 
 


#MBM(F,x,y,fxy,fx1,K)
MBM:=proc(F,x,y,fxy,fx1) local gu,mu,i,lu,gu1,i1:
mu:=subs(y=1,FirstK(F,x,y,fxy,fx1,5)):
mu:=add(mu[i+1]*x^i,i=0..5):
gu:=factor(discrim(discrim(numer(F-fxy),fxy),y)):


for i from 1 to nops(gu) do
 gu1:=op(i,gu):


 if normal(diff(gu1,fx1))<>0 then
  lu:=taylor(subs(fx1=mu,gu1),x=0,7):
   if {seq(expand(coeff(lu,x,i1)),i1=0..5)}={0} then
    gu1:=add(factor(coeff(gu1,fx1,i))*fx1^i,i=0..degree(gu1,fx1)):
    RETURN(gu1):
   fi:
 fi:

od:

print(`The double discriminant is`, gu):
	   
FAIL:


end: