#SCHUTZENBERGER: A PACKAGE FOR HANDLING AND CONJECTURING AND PROVING
#ALGEBRAIC EQUATIONS SATISFIED BY INTEGER SEQUENCES Written by Doron
#Zeilberger, Temple University <zeilberg@math.temple.edu>
 
#LAST UPDATE: Jan. 14, 1997
 
#The most current version is available from 
#http://www.math.temple.edu/~zeilberg
 
#To use it get it first, call it `SCHUTZENBERGER` then go into Maple
#by typing `maple`. Then, in Maple , type `read SCHUTZENBERGER;`
#and see the help files by doing ezra() and ezra(function_name);
 
#(Lots of overlap with Salvy and Zimmerman's Maple package
#gfun available by anon. ftp to inria.fr)
#Contains the following procedures:
 
#difftorec(ope1,D,x,n,N): Converts the differential eq. (linear
#with poly coeff. satisfied by the g.f f(x) to the  recurrence
#satisfied by the coeffs. of f(x)
 
 
#algtodiff(F,P,x,ORDER): Using Comtet's method, given a polynomial
#F(P,x), finds the linear diff. eq. of order ORDER satisfied by the
#formal power series f(x) that satisfies the alg.eq. F(f(x),x)=0 
#Given an algebraic equation F(P(x),x) and desired
#order, finds, if one exists, a linear diff. eq. of order ORDER
 
#genseq(eq,P,x,K) Generates the first K terms of the formal powers series
#that is a solution of the eq P=eq(x,P)
 
 
#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)
 
#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
# gu is a list
 
 
 
 
#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
 
 
#newts(pol,x): gives a list of the power sums of the roots pol(x)
#from p_1 to p_(degree of pol(x)), using newtson's equations
 
 
#newt(pol,x,R): gives a list of the power sums of the roots pol(x)
#from p_1 to p_R, using newtson's equations
 
 
#pol_in_terms_of_p(resh,x) :given the list of power sums [p_1,..,p_deg]
#finds the monic polynomial, in x,  whose roots are the power sums of
 
 
#decomp(pol,x): given a polynomial in x, pol, decides whether it is
#the algebraic product of a quadratic polynomial pol1, and another pol
#pol2. If yes, it returns pol1 and pol2, if not, it returns 0
 
 
#decompx(pol,N,x): given a polynomial in the variables N and x
#investigates whether, when viewed as a polynomial in N
# it can be written as the `alg. product' (i.e. its roots are
#the direct product of the sets of the roots of its `factors')
#of a quadratic in N
#of the form N^2+2*a*x*N-1 where a is an (algebraic) pure number
#and another polynomial in N. with coeffs. that are polynomials in x
 
 print(`SCHUTZENBERGER: A PACKAGE FOR HANDLING AND CONJECTURING AND PROVING`):
 print(`ALGEBRAIC EQUATIONS SATISFIED BY INTEGER SEQUENCES, and RELATED STUFF`):
 print(`Version of Jan. 14, 1997`):
 print(`Named after our great master Marco Schutzenberger(1920-1996)`):
 print(`written by Doron Zeilberger(zeilberg@math.temple.edu).`):
 print(`The most current version is always available from`):
 print(`http://www.math.temple.edu/~zeilberg`):
 print(`(Lots of overlap with Salvy and Zimmerman's Maple package`):
 print(`gfun available by anon. ftp to inria.fr).`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
ezra:=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 `):
 print(`a specific procedure, type ezra(procedure_name)`):
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) finds empirically an ordi. linear recurrence`):
print(` with polynomial coeffs. 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])=`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:
 
 
 
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu:
if (1+DEGREE)*(1+ORDER)+2+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 nops(kv)>1 then
RETURN(0):
fi:
 
for i from 1 to nops(kv) do
  ope:=subs(op(i,kv)=1,ope):
od:
 
ope:
 
end:
 
Findrec:=proc(f,n,N)
local i1,kal,DEGREE, ORDER,L,gu,gu1,i:
 
for L from 1 to nops(f)-5 do
for ORDER from 1 to min(L,trunc((nops(f)-3)/2)-1) do
for DEGREE from 0 to min(L-ORDER,
trunc((nops(f)-3)/2/(1+ORDER))-2) do
gu:=findrec(f,DEGREE,ORDER,n,N):
if gu<>0 then
gu1:=0:
 
for i from 0 to degree(gu,N) do
gu1:=gu1+factor(coeff(gu,N,i))*N^i:
od:
 
kal:=checkrec(gu1,n,N,f):
if kal=[seq(0,i1=1..nops(kal))] then
RETURN(gu1):
fi:
fi:
od:
od:
od:
0:
end:
 
 
#checkrec(ope,n,N,Seq1) applies the operator ope(N,n) to the sequence gu
 
checkrec:=proc(ope,n,N,Seq1)
local i,j,lu,ru:
 
lu:=[]:
 
for i from 1 to nops(Seq1)-degree(ope,N) do
ru:=0:
 
for j from 0 to degree(ope,N) do
ru:=ru+subs(n=i,coeff(ope,N,j))*op(i+j,Seq1):
od:
 
lu:=[op(lu),ru]:
od:
lu:
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,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(0):
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:=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: