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


print(`First Written: March 5, 2024: tested for Maple 2020 `):
print(`Version  : March 5, 2024 `):
print():
print(`This is Tutte.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(` Solving  Functional Equations Dear to  W.T. Tutte using the Naive (yet fullly rigorous!) Guess And Check Method`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/tutte.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Tutte.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():




ezraS:=proc()
if args=NULL then

print(`The Story procedures are:  Info0, InfoA, InfoR, Paper0`):
print(``):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: algtorec, Empir, Findrec, JRratio, SerEx  `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Tutte.txt: A Maple package for  automatically solving functional equations of a type that was needed by W.T. Tutte in his map enumeration studies`):

print(`The MAIN procedures are: GaP, GaPfull, SE, Seq0  `):


print(``):

elif nargs=1 and args[1]=algtorec then
print(`algtorec(F,P,x,n,N): Inputs a polynomial F in the two variables P and x, where P is the ordinary generating function`):
print(`of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator`):
print(`for the Catalan numbers when they are defined by P=1+x*P^2, try:`):
print(`algtorec(P-1-x*P^2,P,x,n,N);`):


elif nargs=1 and args[1]=CheckSerEx then
print(`CheckSerEx(FE,psi,g,x,y,N,R): checks  SerEx(FE,psi,g,x,y,N,R). Try`):
print(`CheckSerEx(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,2);`):


elif nargs=1 and args[1]=Empir then
print(`Empir(L,x,P): guesses an algebraic equation , F(x,P) satisfisfied by P=Sum(L[i]*x^(i-1),i=1..infinity)  but only using the terms of L. Try:`):
print(`Empir([seq(binomial(2*i,i)/(i+1),i=0..30)],x,P);`):



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 coeffs.`):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=GaP then
print(`GaP(FE,psi,g,x,y,N): inputs a polynomial FE in the variables psi,g,x,y, where psi is short for the formal power series psi(x,y), g is short for psi(x,0)`):
print(`x and y are the independent variables, and N is a positive integer for guessing. Tries to discoever and PROVE an algebraic equation for g=g(x). Try:`):
print(`GaP(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30);`):
print(`GaP(y^2*psi^2+(x+x*g*y-2*y-3*y^2)*psi+y-x*g,psi,g,x,y,50);`):


elif nargs=1 and args[1]=GaPfull then
print(`GaPfull(FE,psi,g,x,y,N): inputs a polynomial FE in the variables psi,g,x,y, where psi is short for the formal power series psi(x,y), g is short for psi(x,0)`):
print(`x and y are the independent variables, and N is a positive integer for guessing. Tries to discoever and PROVE an algebraic equation for g=g(x). Also gives the algebraic equation for psi(x,y). Try:`):
print(`GaPfull(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30);`):
print(`GaPfull(y^2*psi^2+(x+x*g*y-2*y-3*y^2)*psi+y-x*g,psi,g,x,y,50);`):


elif nargs=1 and args[1]=Info0 then
print(`Info0(FE,psi,g,x,y,MaxC,n,N,GA): inputs a polynomial FE in the variables psi,g,x,y, representing the functional equation FE(psi,g,x,y)=0 where psi=psi(x,y), g=psi(x,0)`):
print(`and outputs`):
print(`(i) The list of the first few   terms of the coefficients of g(x)=psi(x,0), starting at n=1`):
print(`(ii): a polynomial in x and g=g(x) that gives the algebraic equation satisfied by g(x)`):
print(`(iii) a polynomial in x,y, and psi=psi(x,y) giving the algebraic equation satisfied by psi(x,y)`):
print(`(iv) the pair [INI,ope] where INI is the list of initial conditions and ope is the linear recurrence operator with constant coefficients`):
print(`annihilationg the sequence of coefficients of g(x). For W.T. Tutte's functional equation in his seminal paper, type:`):
print(`(v): The coeff. of x^GA in g(x) `):
print(`Info0(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,10,n,N,1000);`):

elif nargs=1 and args[1]=InfoA then
print(`InfoA(FE,psi,g,x,y,N,R,P), guessed algeraic equations for P[i] where P[i] is the  coeff. of  y^(i-1) in psi=psi(x,y), where psi(x,y) solution of the`):
print(`functional equation FE(psi,g)=- where g=psi(x,0)`):
print(`InfoA(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,3,P);`):

elif nargs=1 and args[1]=JRratio then
print(`JRratio(FE,psi,g,x,y,N,r): The sequence of ratios of the coefficient of x^n in (g(x)-1)^r divided by the coeff. of x^n of g(x). Try:`):
print(`JRratio(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,2);`):


elif nargs=1 and args[1]=Paper0 then
print(`Paper0(FE,psi,g,x,y,MaxC,n,GA): inputs a polynomial FE in the variables psi,g,x,y, representing the functional equation FE(psi,g,x,y)=0 where psi=psi(x,y), g=psi(x,0)`):
print(`and outputs a paper with`):
print(`(i) The list of the first few   terms of the coefficients of g(x)=psi(x,0), starting at n=1`):
print(`(ii): a polynomial in x and g=g(x) that gives the algebraic equation satisfied by g(x)`):
print(`(iii) a polynomial in x,y, and psi=psi(x,y) giving the algebraic equation satisfied by psi(x,y)`):
print(`(iv) the pair [INI,ope] where INI is the list of initial conditions and ope is the linear recurrence operator with constant coefficients`):
print(`annihilationg the sequence of coefficients of g(x). For W.T. Tutte's functional equation in his seminal paper, type:`):
print(`(v): The coeff. of x^GA in g(x) `):
print(`Paper0(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,5,n,1000);`):

elif nargs=1 and args[1]=SE then
print(``):
print(`SE(FE,psi,g,x,y,K): the first K+1 terms of of the series expansion, in x,  of the solution of the functional equation`):
print(`FE=0, where FE is a polynomial in psi=psi(x,y) and g=g(x), where psi denotes a formal power series in (x,y), psi(x,y)`):
print(`and g is short for psi(x,0). Try:`):
print(`SE(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30);`):

elif nargs=1 and args[1]=Seq0 then
print(`Seq0(FE,psi,g,x,y,N,R): the list of length R+1 whose i-th entry is the list of the first N coefficients, of the coeff. of y^(i-1)starting at n=1 of SE(FE,psi,g,x,y,N). Try:`):
print(`Seq0(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,3);`):

elif nargs=1 and args[1]=SerEx then
print(`SerEx(FE,psi,g,x,y,N,R): Outputs GaP(FE,psi,g,x,y,N) followed by the expressions in g of the coeffs. of y^i SE(FE,psi,g,x,y,N) up to i=R`. Try):
print(`SerEx(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,2);`):

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

fi:


end:




#Start From FindRec.txt



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,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 (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:

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

#End From FindRec.txt

###START FROM SCHUTZENBERGER.txt
 

#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:
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(FAIL):
fi:
 
flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=F/pip:
normal(F):
end:
 
#Empir(L,x,P): guesses an algebraic equation , F(x,P) satisfisfied by P=Sum(L[i]*x^(i-1),i=1..infinity)  but only using the terms of L. Try:
#Empir([seq(2*i.i)/(i+1),i=0..30)],x,P);

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<>FAIL then
RETURN(lu):
fi:
od:
od:
od:
FAIL:
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:
 
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:

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

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

ope:=subs(var,ope):

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

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:

 
pashet:=proc(p,n,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:

#MyGcd(L): inputs a list of polynomials and finds their gcd. Try
#MyGcd([n,n^2,n^3]);
MyGcd:=proc(L) local gu:
if L=[] then
 RETURN(FAIL):
elif nops(L)=1 then
 RETURN(L[1]):
else
 gu:=MyGcd([op(1..nops(L)-1,L)]):
 RETURN(gcd(gu,L[nops(L)])):
fi:
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,N):
end:


#algtorec(F,P,x,n,N): Inputs a polynomial F in the two variables P and x, where P is the ordinary generating function
#of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator
#for the Catalan numbers when they are defined by P=1+x*P^2, try
#algtorec(P-1-x*P^2,P,x,n,N);
algtorec:=proc(F,P,x,n,N) local ope,d,i, mu:

for d from 1 to degree(F,P) do
ope:=algtodiff(F,P,x,d):
 if ope<>FAIL then
   ope:=difftorec(ope,D,x,n,N):
  mu:=MyGcd([seq(coeff(ope,N,i),i=0..degree(ope,N))]):
  ope:=normal(ope/mu):
  ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
 RETURN(ope):   
 fi:
od:
FAIL:
end:

 
###END FROM SCHUTZENBERGER.txt


#SE(FE,psi,g,x,y,N): the first N+1 terms of of the series expansion of the solution of the functional equation
#FE=0, where FE is a polynomial in psi and g and psi denotes a formal power series in (x,y), psi(x,y)
#and g is short for psi(x,0). Try
#SE(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30);
SE:=proc(FE,psi,g,x,y,N) local f,kha,kha0,kha1,lu,n:
option remember:
f:=1:

for n from 1 to N do
if subs(y=0,denom(f))=0 then
 RETURN(FAIL):
fi:
lu:=expand(subs({psi=f+x^n*kha,g=subs(y=0,f+x^n*kha0)},FE));
lu:=coeff(lu,x,n):
if normal(diff(lu,kha0))<> 0 then
 RETURN(FAIL):
fi:

kha1:=normal(solve(lu,kha)):
f:=f+x^n*kha1:
od:

f:
end:



#GaP(FE,psi,g,x,y,N): inputs a polynomial FE in the variables psi,g,x,y, where psi is short for the formal power series psi(x,y), g is short for psi(x,0)
#x and y are the independent variables, and N is a positive integer for guessing. Tries to discoever and PROVE an algebraic equation for g=g(x). Try:
#GaP(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30);
GaP:=proc(FE,psi,g,x,y,N) local lu,eq,ku,i:

lu:=SE(FE,psi,g,x,y,N): 
if lu=FAIL then
 RETURN(FAIL):
fi:

lu:=subs(y=0,lu):

lu:=[seq(coeff(lu,x,i),i=0..N)]:

eq:=Empir(lu,x,g):

if eq=FAIL then
 print(` Make `, N, `bigger `):
 RETURN(FAIL):
fi:

ku:=eliminate({eq,FE},g)[2][1]:

lu:=[eq,subs(psi=g,factor(subs(y=0,ku)))]:


if diff(denom(normal(lu[2]/lu[1])),g)<>0 then
print(`The conjectured equation`):
print(eq):
print(`didn't work out`):

 RETURN(FAIL):
fi:

collect(eq,g):

end:





#GaPfull(FE,psi,g,x,y,N): inputs a polynomial FE in the variables psi,g,x,y, where psi is short for the formal power series psi(x,y), g is short for psi(x,0)
#x and y are the independent variables, and N is a positive integer for guessing. Tries to discoever and PROVE an algebraic equation for g=g(x). Try:
#`GaPfull(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30);
GaPfull:=proc(FE,psi,g,x,y,N) local lu,eq,ku,i,ku1,ku2:

lu:=SE(FE,psi,g,x,y,N): 
if lu=FAIL then
 RETURN(FAIL):
fi:

lu:=subs(y=0,lu):

lu:=[seq(coeff(lu,x,i),i=0..N)]:

eq:=Empir(lu,x,g):


if eq=FAIL then
 print(` Make `, N, `bigger `):
 RETURN(FAIL):
fi:


ku:=eliminate({eq,FE},g)[2][1]:


ku1:=eliminate({ku,FE},psi)[2][1]:

ku2:=subs(psi=g,subs(y=0,ku)):

if normal(diff(ku1/eq,g))<>0 then
 RETURN(FAIL):
fi:


if normal(diff(ku2/eq,g))<>0 then
 RETURN(FAIL):
fi:


[eq,ku]:
end:





#JRratio(FE,psi,g,x,y,N,r): The sequence of ratios of the coefficient of x^n in (g(x)-1)^r divided by the coeff. of x^n of g(x). Try:
#JRratio(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,2);
JRratio:=proc(FE,psi,g,x,y,N,r) local L,f,f1,i:

f:=subs(y=0,SE(FE,psi,g,x,y,N)):
f:=f-subs(x=0,f):

f1:=f^r:
L:=[seq(coeff(f1,x,i)/coeff(f,x,i),i=1..N)]:
L:
end:


#SerEx(FE,psi,g,x,y,N,R): Outputs GaP(FE,psi,g,x,y,N) followed by the expressions in g of the coeffs. of y^i SE(FE,psi,g,x,y,N) up to i=R
SerEx:=proc(FE,psi,g,x,y,N,R)  local gu,ku,lu,g1,G,lu1,i,i1:
gu:=GaPfull(FE,psi,g,x,y,N):
ku:=gu[1]:
lu:=gu[2]:
G[0]:=g:
for i from 1 to R do
lu1:=coeff(expand(subs(psi=add(G[i1]*y^i1,i1=0..i-1)+g1*y^i,lu)),y,i);


G[i]:=normal(solve(lu1,g1)):

G[i]:=rem(numer(G[i]),ku,g)/rem(denom(G[i]),ku,g):
od:

[seq(G[i],i=1..R)]:
end:

#CheckSerEx(FE,psi,g,x,y,N,R): checks SerEx(FE,psi,g,x,y,N,R).
CheckSerEx:=proc(FE,psi,g,x,y,N,R) local gu,lu,i,g0,fu,i1:
lu:=SerEx(FE,psi,g,x,y,N,R):
gu:=SE(FE,psi,g,x,y,N):
gu:=taylor(gu,y=0,R+1):
g0:=subs(y=0,gu):

for i from 1 to R do
fu:=taylor(subs(g=g0,lu[i])-coeff(gu,y,i),x=0,N+1):
if {seq(coeff(fu,x,i1),i1=0..N-4)}<>{0} then
 RETURN(false):
fi:
od:
true:
end:



#Seq0(FE,psi,g,x,y,N): the list of length R+1, 
#whose r-th entry is the list of the first N coefficients, starting at n=1 of  to n=N, of the coeffient of y^(r-1),SE(FE,psi,g,x,y,N). Try:
#Seq0(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,3);

Seq0:=proc(FE,psi,g,x,y,N,R) local f,i,r,gu:
gu:=SE(FE,psi,g,x,y,N):
if gu=FAIL then
  RETURN(FAIL):
fi:
f:=taylor(gu,y=0,R+1):


[seq(
[seq(coeff(normal(coeff(f,y,r)),x,i),i=1..N)],
r=0..R)]:

end:



#InfoA(FE,psi,g,x,y,N,R,P), guessed algeraic equations for P[i] where P[i] is the  coeff. of  y^(i-1) in psi=psi(x,y), where psi(x,y) solution of the
#functional equation FE(psi,g)=- where g=psi(x,0)
#InfoA(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,3,P);
InfoA:=proc(FE,psi,g,x,y,N,R,P) local f,lu,i,lu1,f1,L1,i1:

f:=SE(FE,psi,g,x,y,N):

f:=taylor(f,y=0,R+1):

lu:=[]:

for i from 0  to R do
f1:=normal(coeff(f,y,i)):
f1:=taylor(f1,x=0,N+1):
L1:=[seq(coeff(f1,x,i1),i1=0..N)]:

 lu1:=Empir(L1,x,P[i]):
 if lu1=FAIL then
   RETURN(lu):
 else
 lu:=[op(lu),lu1]:
 fi:
od:

lu:

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:
 

#Info0(FE,psi,g,x,y,MaxC,n,N,GA): inputs a polynomial FE in the variables psi,g,x,y, representing the functional equation FE(psi,g,x,y)=0 where psi=psi(x,y), g=psi(x,0)
#and outputs
#(i) The list of the first few terms of the coefficients of g(x)=psi(x,0)
#(ii): a polynomial in x and g=g(x) that gives the algebraic equation satisfied by g(x)
#(iii) a polynomial in x,y, and psi=psi(x,y) giving the algebraic equation satisfied by psi(x,y)
#(iv) the pair [INI,ope] where INI is the list of initial conditions and ope is the linear recurrence operator with constant coefficients
#(v): the coeff. of x^GA in g(x)=psi(x,0)
#annihilationg the sequence of coefficients of g(x). For W.T. Tutte's functional equation in his seminal paper, type:
#Info0(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,n,N);
Info0:=proc(FE,psi,g,x,y,MaxC,n,N,GA) local L,ku,gu,x1,ope,K,Ini:

K:=max(seq((2+MaxC-x1)*(1+x1)+4,x1=0..MaxC))+4:


L:=Seq0(FE,psi,g,x,y,K,0):

if L=FAIL  then
 RETURN(FAIL):
fi:
L:=L[1]:

gu:=[L]:


ku:=GaPfull(FE,psi,g,x,y,K):

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

gu:=[L, op(ku)]:

ope:=Findrec(L,n,N,MaxC):

if ope=FAIL then
 RETURN(gu):
else
Ini:=[op(1..degree(ope,N),L)];
 gu:=[op(gu),[Ini,ope]]:
fi:

[op(gu), SeqFromRec(ope,n,N,Ini,GA)[GA]]:

end:



#Paper0(FE,psi,g,x,y,MaxC,n,GA): inputs a polynomial FE in the variables psi,g,x,y, representing the functional equation FE(psi,g,x,y)=0 where psi=psi(x,y), g=psi(x,0)
#and outputs a paper about
#(i) The list of the first few terms of the coefficients of g(x)=psi(x,0)
#(ii): a polynomial in x and g=g(x) that gives the algebraic equation satisfied by g(x)
#(iii) a polynomial in x,y, and psi=psi(x,y) giving the algebraic equation satisfied by psi(x,y)
#(iv) the pair [INI,ope] where INI is the list of initial conditions and ope is the linear recurrence operator with constant coefficients
#(v): the coeff. of x^GA in g(x)=psi(x,0)
#annihilationg the sequence of coefficients of g(x). For W.T. Tutte's functional equation in his seminal paper, type:
#Paper0(y^2*psi^2+(x+x*g*y-y-y^2)*psi+y-x*g,psi,g,x,y,30,n);

Paper0:=proc(FE,psi,g,x,y,MaxC,n,GA) local N,gu, Psi,a,i,Ini,ope,t0:
t0:=time():
gu:=Info0(FE,psi,g,x,y,MaxC,n,N,GA):

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

print(`On the coefficients of Psi(x,0) of the unique solution of the functional equation`, subs({psi=Psi(x,y),g=Psi(x,0)},FE)=0 ):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let's g(x)=Psi(x,0). The first`, nops(gu[1]), `(for the sake of the OEIS) are `):
print(``):
lprint(gu[1]):
print(``):
if nops(gu)>=2 then
print(`g=g(x) satisfies the following algebraic equation `):
print(``):
print(gu[2]=0):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu[2]=0):
print(``):
print(`More generally, psi=Psi(x,y) satisisfies the algebraic equation `):
print(``):
print(gu[3]=0):
print(``):
print(`and in Maple  notation`):
print(``):
lprint(gu[3]=0):
print(``):
fi:

if nops(gu)>=4 then
 Ini:=gu[4][1]:
 ope:=gu[4][2]:
print(`writing Psi(x,0)=g(x) as a Taylor series around x=0`):
print(``):
print(g(x)= Sum(a[n]*x^n,n=0..infinity)):
print(``):
print(`The coefficients, a[n], satisfy the folllowing linear recurrence equation with polynomial coefficients of order`, nops(Ini)):
print(``):
print(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N))=0):
print(``):
print(`subject to the initial conditions `):
print(``):
lprint(seq(a(i)=Ini[i],i=1..nops(Ini))):
print(``):
print(`Finally, just for fun here is `, a(GA) ):
print(``):
lprint( a(GA)=gu[5]):
print(``):
fi:
print(``):
print(`---------------------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `seconds to generate `):
print(``):

end: