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


print(`First Written: Feb. 2022: tested for Maple 2020 `):
print():
print(`This is PureRecRat.txt, one of  Maple packages`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`Linear-Time and Constant-Space Algorithms to compute Multi-Sequences that arise in Enumerative Combinatorics  `):

with(combinat):

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

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

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

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

print(`For a list of the STORY functions type: ezraS();`):
print(` For specific help type "ezra(procedure_name);" `):
print():

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




ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` CheckScheme, CoefD, CoefDc, GJgf, INIc, INItay,  SchemeD, Walk2Ddirect `):
else
ezra(args):

fi:


end:


ezraS:=proc()
if args=NULL then

print(`The STORY procedures are`):
print(` BookWalks,SchemeV, Walk2Dpaper`):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` PureRec.txt: A Maple package for computing with pure recurrence relations  using the Apagodu-Zeilberger multivariate extension of the Almkvist-Zeilberger algorithm `):

print(`The MAIN procedures are`):
print(` EvalScheme, Scheme, SchemeC `):


elif nargs=1 and args[1]=BookWalks then
print(`BookWalks(N1,n1,K):  An article with lots of theorerms about lattice walks with sets of atomic stape of size from 2 to n1, and coordinates up to N1. `):
print(`K is as in Walk2Dpaper(St,K) (q.v.) . Try:`):
print(`BookWalks(2,3,1000):`):

elif nargs=1 and args[1]=CheckScheme then
print(`CheckScheme(R,n,x,K1,K2): Checks Scheme(R,n,x,m,M) for K2 different values in [0,K1]^n. Try:`):
print(`CheckScheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,10,10);`):

elif nargs=1 and args[1]=CoefD then
print(`CoefD(F,x,m): Inputs a function of x[1], ..., x[n] (where m:=nops(m)) that has a nice Taylor expansion around the origin, and a list m of n (say) non-neg. integers, outputs`):
print(`the Taylor coefficient (around the origin, in other words, the MacLaurin coefficient) of x[1]^m[1]*...*x[n]^m[n]. Try:`):
print(`CoefD(1/(1-x[1]-x[2]-x[3]-x[4]+x[1]*x[2]*x[3]*x[4]),x,[2,3,4,5]);`):




elif nargs=1 and args[1]=EvalScheme then
print(`EvalScheme(S,m,M,m0): inputs a scheme S (let n:=nops(S)) in n dimension, phrased in terms of m and M (i.e. the discrete variables are m[1],...,m[n] and the corresponding shift operators are`):
print(`M[1],..., M[n], and a list of non-negative integers m0. Evaluats it at m0. Try:`):

print(`EvalScheme([ [2/M[1],3/M[2],4/M[3]],[[[1]]],m,M,[2,2,2]);`):

print(`EvalScheme([[1/M[1]+1/M[1]^2],[1,1]],m,M,[3]);`):



elif nargs=1 and args[1]=GJgf then
print(`GJgf(alphabet,MISTAKES,x,t): The Goulden-Jackson generating function`):
print(` GJgf({1,2},{[1,2,1],[2,1,2]},x,t); `):

elif nargs=1 and args[1]=GJpaper then
print(`GJpaper(MIS,K): inputs a set of words, MIS, in {1,2} and outputs a paper about the number of words in 1^a2^b avoiding the members of M as consectuive wubsord.. Try:`):
print(`Still needs work. Do NOT use`):
print(`GJpaper({[1,2,1],[2,1,2]},1000);`):

elif nargs=1 and args[1]=INIc then
print(`INIc(P,x,L): Given a Laurent polynomial in the list of variables x, (of size k, say) and a list L of size k+1, outputs`):
print(`a list-of-lists-....L[m1][m2]..[m[k+1]] is the coefficient of x[1]^m[1]*...*x[k]^m[k] in P^m[k+1]`):
print(`The initial conditions up. Try: `):
print(`INIc(x[1]+1+x[1],[x[1]],[2,3]);`);

elif nargs=1 and args[1]=INItay then
print(`INItay(R,x,A): Given a function R in the list of variables x[1],..., x[n] (for some n:=nops(x)) finds the list of lists...of lists, let's call it L, such that for each 0<=a[1]<A[1], ..., 0<=a[n<A[n]`):
print(`The coefficient of x[1]^a[1]*...*x[n]^a[n] in the Taylor expansion is L[a[1]+1]...[a[n]+1]. Try:`):
print(`INItay(1/(1-x-y-z+x^2*y^2+z^2),[x,y,z],[3,3,3]);`):


elif nargs=1 and args[1]=Scheme then
print(`Scheme(R,k,x,m,M): Let k be a positive integer (the second input parameter) and let R be a function of x[1],..., x[k]`):
print(`that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of`):
print(`the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m[1],...,m[k] `):
print(`and the respective shift operators M[1],...,M[k]. Try:`):
print(`Scheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,m,M);`):


elif nargs=1 and args[1]=SchemeC then
print(`SchemeC(P,k,x,m,M,n,N): Let k be a positive integer (the second input parameter) and let P be a Laurent polynomial in x[1],..., x[k]`):
print(`uses the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m[1],...,m[k], m[k+1]`):
print(`and the respective shift operators M[1],...,M[k], M[k+1], for the discrete function F(m[1],..., m[k+1]):=`):
print(`The coeff. of x[1]^m[1]*...*x[k]^m[k] in P^m[k+1]. Try:`):
print(`SchemeC(x[1]+1/x[1]+1+x[2]+1/x[2],2,x,m,M);`):



elif nargs=1 and args[1]=SchemeD then
print(`SchemeD(R,n,x,m,M): Let n be a positive integer (the second input parameter) and let R be a function of x[1],..., x[n]`):
print(`that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of`):
print(`the Almkvist-Zeilberger algorithm to output a PURE recurrence in terms of the discrete variables m`):
print(`and the  shift operator M for the coefficient of (x[1]*...*x[n])^m`):
print(`SchemeD(1/(1-x[1]-x[2]-x[3]-4*x[1]*x[2]*x[3]),3,x,m,M);`):

elif nargs=1 and args[1]=SchemeV then
print(`SchemeV(R,n,x): A verbose version of Scheme(R,n,x,m,M) with illustration. try:`):
print(`SchemeV(1/(1-x[1]-x[2]-x[1]*x[2])^(1/3),2,x);`):


elif nargs=1 and args[1]=Walk2Ddirect then
print(`Walk2Ddirect(St,pt): The number of ways of walking from [0,0] to pt using the set of steps St, done directly. Only for checking. Try:`):
print(`Walk2Ddirect({[1,0],[0,1],[1,1]},[20,30]);`):

elif nargs=1 and args[1]=Walk2Dpaper then
print(`Walk2Dpaper(St,K): inputs a set of positive steps,St and outputs an article about the bi-sequence: Number of walks from the origin to a point [m1,m2].`):
print(`It gives the pure recurrence scheme  and computes as an illustrative example`):
print(`the number of ways of walking to [2*K,K] and [K,K]. It also gives the recurrence for the diagonal (i.e. the one-variable  and recurrence) for getting from the origin to [0,0] to [n,n]. For example, for`):
print(`the Royal walks : try:`):
print(`Walk2Dpaper({[1,0],[0,1],[1,1]},1000);`):

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

fi:


end:




#FROM David_Ian 


 
#findeq sets up the equ C[v]= x[v]+t*Sum_u overlap(u,v,x) *C[u]




hakten:=proc(B)
local w,i:
 
for i from 1 to nops(B) do
w:=op(i,B):
 
if superflous(B,w)=1 then
 RETURN(B minus {w}):
fi:
od:
B:
end:

#Hakten(B) removes all superflous words
 
Hakten:=proc(B)
local B1,B2:
 
B1:=B:
 
B2:=hakten(B1):
 
while B1<>B2 do
B1:=B2:
B2:=hakten(B2):
od:
 
B2:
end:

#overlap is a procedure that given two words u and v
#computes the weight-enumerator of all v\suffix(u),
#for all suffixes of u that are prefixes of v.
overlap:=proc(u,v,x)
local i,j,k,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
  for k from j-i+1 to nops(v) do
   gug:=gug*x[op(k,v)]:
  od:
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:
 
findeq:=proc(v,MISTAKES,C,t,x)
local eq,i,u:
eq:=t:
 
for i from 1 to nops(v) do
 eq:=eq*x[op(i,v)]:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq+t*overlap(u,v,x)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:

GJgf:=proc(alphabet,MISTAKES,x,t)
local v,eq, var,i,lu,C:
 
if Hakten(MISTAKES)<>MISTAKES then
ERROR(`The set of mistakes is not minimal, use another package`):
fi:
 
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeq(v,MISTAKES,C,t,x)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 
for i from 1 to nops(alphabet) do
 lu:=lu-x[op(i,alphabet)]:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
lu:=normal(1/lu):
subs(t=t-1,lu): 
end:
 

#End FROM David_Ian 

#START FROM MultiAlmkvistZeilberger

 
ezraMAZ:=proc()
if nargs=0 then
print(`MultiAlmkvistZeilberger`):
print(`A Maple package for  finding recurrences for multi-integrals of `):
print(` Hypergeometric/Hyperexponential Integrands.`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  MAZ, MAZpaper `):
 
elif nargs=1 and args[1]=MAZ then
print(`MAZ(POL,H,rat,x,n,N,pars): Given a polynomial, POL, an hyper-exponential function, H, and`):
print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`):
print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`):
print(`ouputs a recurrence operator, let's call it  ope(n,N), and a multi-certificate, let's call it R`):
print(`such that the function F(n;x):=POL*H*rat^n satisfies`):
print(` ope(n,N)F=sum(D_{x[i]} (R[i]H*rat^n),i=1..nops(R)`):
print(`In particular, try:`):
print(`MAZ(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,A,{b,c});`):
print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{});`):



elif nargs=1 and args[1]=MAZpaper then
print(`MAZpaper(POL,H,rat,x,n,pars,A,R): Given a polynomial, POL, an hyper-exponential function, H, and`):
print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`):
print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`):
print(`and letters A and R for expressing the statement of the theorem`):
print(`and proof`):
print(`ouputs a paper stating and proving `): 
print(`A recurrence equation for the (contour or vanishing-at-endpoints) integral`):
print(`of F:=POL*H*rat^n `):
print(`with respect to the variables in the list [x] `):
print(`Try:`):
print(`MAZpaper(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,{b,c},A,R);`):
print(`MAZpaper(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,{},A,R);`):


else print(`There is no help for`, args):

fi:



end:



#FROM David_Ian 


 
#findeq sets up the equ C[v]= x[v]+t*Sum_u overlap(u,v,x) *C[u]

ez:=proc(): print(` GJgf(alphabet,MISTAKES,x,t) `):end:



hakten:=proc(B)
local w,i:
 
for i from 1 to nops(B) do
w:=op(i,B):
 
if superflous(B,w)=1 then
 RETURN(B minus {w}):
fi:
od:
B:
end:

#Hakten(B) removes all superflous words
 
Hakten:=proc(B)
local B1,B2:
 
B1:=B:
 
B2:=hakten(B1):
 
while B1<>B2 do
B1:=B2:
B2:=hakten(B2):
od:
 
B2:
end:

#overlap is a procedure that given two words u and v
#computes the weight-enumerator of all v\suffix(u),
#for all suffixes of u that are prefixes of v.
overlap:=proc(u,v,x)
local i,j,k,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
  for k from j-i+1 to nops(v) do
   gug:=gug*x[op(k,v)]:
  od:
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:
 
findeq:=proc(v,MISTAKES,C,t,x)
local eq,i,u:
eq:=t:
 
for i from 1 to nops(v) do
 eq:=eq*x[op(i,v)]:
od:
 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq+t*overlap(u,v,x)*C[op(u)]:
od:
 
C[op(v)]-eq=0:
 
end:

GJgf:=proc(alphabet,MISTAKES,x,t)
local v,eq, var,i,lu,C:
 
if Hakten(MISTAKES)<>MISTAKES then
ERROR(`The set of mistakes is not minimal, use another package`):
fi:
 
eq:={}:
var:={}:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 eq:= eq union {findeq(v,MISTAKES,C,t,x)}:
 var:=var union {C[op(v)]}:
od:
 
var:=solve(eq,var):
 
lu:=1:
 
for i from 1 to nops(alphabet) do
 lu:=lu-x[op(i,alphabet)]:
od:
 
for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 lu:=lu-subs(var,C[op(v)]):
od:
 
lu:=normal(1/lu):
subs(t=t-1,lu): 
end:
 

#End FROM David_Ian 



_EnvAllSolutions:=true:

####From MultiZeilberger
#IV1(d,n): all the vectors of non-negative integres of length d whose
#sum is n
IV1:=proc(d,n)
local gu,i,gu1,i1:
if d=0 then
 if n=0 then
   RETURN({[]}):
 else
    RETURN({}):
fi:
fi:
 
gu:={}:

for i from 0 to n do
 gu1:=IV1(d-1,n-i):
 gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}:
od:
gu:
end:

yafe:=proc(ope,N) local i:
add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):end:

#IV(d,n): all the vectors of non-negative integres of length d whose
#sum is <=n
IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end:

#GenPol(kList,a,deg): The generic polynomial in kList of
#degree deg,  using the indexed variable a, followed by the set
#of coeffs.
GenPol:=proc(kList,a,deg)
local gu,i,i1:
gu:=IV(nops(kList),deg):
convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`),
i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}:

end:

#Extract1(POL,kList): extracts all the coeffs. of a POL in kList
Extract1:=proc(POL,kList) local k1,kList1,i:
k1:=kList[nops(kList)]:
kList1:=[op(1..nops(kList)-1,kList)]:
if nops(kList)=1 then
 RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}):

else
 RETURN({seq(
op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}):
fi:

end:

###########End from MultiZeilberger


FindDeg:=proc(POL,H,rat,x,xSet,n,L)
local s,t,h,e,i,Hbar,q,r,gu:
s:=numer(rat): t:=denom(rat):

h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`):

Hbar:=H*s^n/t^(n+L):

gu:=normal(simplify(diff(Hbar,x)/Hbar)):
q:=numer(gu): r:=denom(gu):
max(degree(h,xSet)-degree( diff(r,x)+q,xSet)+degree(r,xSet) ,degree(h,xSet)-degree(r,xSet)+1+degree(r,xSet)):
end:

#MAZ1(POL,H,rat,x,n,N,L,pars): Given a polynomial POL of the discrete variable n 
#and the continuous variables in x
#and a pure-hyperexponential function H of the variables in x, and a rational function
#Rat of x, and a shift operator N, and a non-neg. integer L
#returns FAIL or, if there is one,  an operator ope(N,n) of order L and
#a list of rational functions R, such that F:=POL*H*rat^n satisfies
#ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) .
#For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{})
MAZ1:=proc(POL,H,rat,x,n,N,L,pars)
local deg,deg1,m,gu,i,i1:
deg:=[]:
for i from 1 to nops(x) do
deg:=[op(deg),FindDeg(POL,H,rat,x[i],{seq(x[i1],i1=1..nops(x))},n,L)]:
od:

m:=min(op(deg)):

for i from 0 to m do
 deg1:=[seq(deg[i1]-m+i,i1=1..nops(deg))]:

 gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,deg1):

  if gu<>FAIL then
     RETURN(gu):
  fi:
od:
FAIL:
end:


MAZ:=proc(POL,H,rat,x,n,N,pars) local gu,L:

for L from 0 do
gu:=MAZ1(POL,H,rat,x,n,N,L,pars):

if gu<>FAIL then
RETURN(gu):
fi:
od:
end:


CheckMAZ:=proc(POL,H,rat,x,n,N,pars)
local F,gu,luL,luR,R,ope,i:
F:=H*rat^n:
gu:=MAZ(POL,H,rat,x,n,N,pars):
if gu=FAIL then
RETURN(FAIL):
fi:
ope:=gu[1]:
R:=gu[2]:
luL:=normal(add(subs(n=n+i,POL)*coeff(ope,N,i)*normal(simplify(subs(n=n+i,F)/F)),i=0..degree(ope,N))):
luR:=normal(add( normal(diff(R[i]*F,x[i])/F ),
              i=1..nops(R))):
evalb(normal(luL/luR)=1):
end:


#Walk2Ddirect(St,pt): The number of ways of walking from [0,0] to pt using the set of steps St, done directly. Only for checking. Try:
#walk2Ddirect({[1,0],[0,1],[1,1]},20,30);
Walk2Ddirect:=proc(St,pt) local s:
option remember:

if pt[1]<0 or pt[2]<0 then
0:
elif pt=[0,0] then
1:
else
add(Walk2Ddirect(St,pt-s) , s in St):
fi:
end:



#MAZ1deg(POL,H,rat,x,n,N,L,pars,deg): Given a polynomial POL of the discrete variable n 
#and the continuous variables in x
#and a pure-hyperexponential function H of the variables in x, and a rational function
#Rat of x, and a shift operator N, and a non-neg. integer L
#returns FAIL or, if there is one,  an operator ope(N,n) of order L and
#a list of rational functions R, such that F:=POL*H*rat^n satisfies
#ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) .
#For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{})
MAZ1deg:=proc(POL,H,rat,x,n,N,L,pars,deg)
local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,ku,i1,eqN,var1N,opeN,bilti,meka:

s:=numer(rat): t:=denom(rat):
ope:=add(e[i]*N^i,i=0..L):
h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`):

Hbar:=H*s^n/t^(n+L):
gu:=[]:

for i from 1 to nops(x) do
 gu:=[op(gu),normal(simplify(diff(Hbar,x[i])/Hbar))]:
od:

q:=[]: r:=[]:
for i from 1 to nops(gu) do
q:=[op(q),numer(gu[i])]: 
r:=[op(r),denom(gu[i])]:
od:

var:={}:
X:=[]:
for i from 1 to nops(x) do
 ku:=GenPol(x,a[i],deg[i]):
 X:=[op(X),ku[1]]:
 var:=var union ku[2]:
od:

 var:=var union {seq(e[i],i=0..L)}:

gu:=h:

for i from 1 to nops(x) do
gu:=expand(normal(gu-(diff(r[i],x[i])+q[i])*X[i]/r[i]-r[i]*diff(X[i]/r[i],x[i]))):
od:
eq:=Extract1(numer(gu),x):

eqN:=subs(n=9/17,eq):
eqN:=subs({seq(pars[i]=1/(i^2+3),i=1..nops(pars))},eqN):

var1N:=solve(eqN,var):

opeN:=subs(var1N,ope):

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



var1:=solve(eq,var):

ope:=subs(var1,ope):
X:=subs(var1,X):

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


bilti:={}:

for i from 1 to nops(var1)  do

 if op(1,op(i,var1))=op(2,op(i,var1)) then
   bilti:=bilti union {op(1,op(i,var1))}:
 fi:
  
od:




for i from 1 to nops(bilti) do
 if coeff(ope,bilti[i],1)<>0 then
    ope:=coeff(ope,bilti[i],1):
    X:=[seq(coeff(X[i1],bilti[i],1),i1=1..nops(X))]:
   
    meka:=coeff(ope,N,degree(ope,N)):
     ope:=expand(normal(ope/meka)):
     ope:=yafe(ope,N):
     X:=[seq(normal(X[i1]/meka/t^L),i1=1..nops(X))]: 
    RETURN(ope,X):
 fi:
od:

end:



#MAZ(POL,H,rat,x,n,pars,F,A,R)
MAZpaper:=proc(POL,H,rat,x,n,pars,A,R) local gu,F,F1,ope,i,N:

print(``):

F1:=H*rat^n:

gu:=MAZ(POL,H,rat,x,n,N,pars):
ope:=gu[1]:
if gu=FAIL then
 RETURN(FAIL):
fi:
print(`A Linear Recurrence Equation for an Integral in`, nops(x), ` variables `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let `):
print(F(n,op(x))=F1):
print(`and let`, A(n), `be its definite integral  of`):
print(POL*F(n,op(x))):
print(`with respect to the continuous`):
print(` variables `, op(x)):
print(``):
print(A(n) , `satisfies the following homogeneous linear recurrence with`):
print(`polynomial coefficients: `):
print(``):
print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Proof: We cleverly constuct the rational functions`):

for i from 1 to nops(gu[2]) do
print(R[i](n,op(x))=gu[2][i]):
od:
print(`With the motive that`):
print(add(coeff(ope,N,i)*F(n+i,op(x)),i=0..degree(ope,N))=
add(D[x[i]](R[i](n,op(x))*F(n,op(x))),i=1..nops(x))):
print(``):
print(`(Check!)`):
print(`and the theorem follows upon integrating with respect to`, op(x)):
print(` QED! `):



end:

###END FROM MultiAlmkvistZeilberger




#CoefD(F,x,m): Inputs a function of x[1], ..., x[n] (where m:=nops(m)) that has a nice Taylor expansion around the origin, and a list m of n (say) non-neg. integers, outputs
#the Taylor coefficient (around the origin, in other words, the MacLaurin coefficient) of x[1]^m[1]*...*x[n]^m[n]. Try:
#CoefD(1/(1-x[1]-x[2]-x[3]-x[4]+x[1]*x[2]*x[3]*x[4]),x,[2,3,4,5]);

CoefD:=proc(F,x,m) local n,i,F1:

n:=nops(m):


F1:=F:
for i from 1 to n do
F1:=taylor(F1,x[i]=0,m[i]+2):
F1:=expand(normal(coeff(F1,x[i],m[i]))):
od:
F1:

end:




#INItay(R,x,A): Given a function R in the list of variables x[1],..., x[n] (for some n:=nops(x)) finds the list of lists...of lists, let's call it L, such that for each 0<=a[1]<A[1], ..., 0<=a[n<A[n]
#The coefficient of x[1]^a[1]*...*x[n]^a[n] in the Taylor expansion is L[a[1]+1]...[a[n]+1]. Try:
#INItay(1/(1-x-y-z+x^2*y^2+z^2),[x,y,z],[3,3,3]);
INItay:=proc(R,x,A) local n,i,R1:
n:=nops(x):
if not (type(x,list) and nops(A)=n) then
 RETURN(FAIL):
fi:

R1:=taylor(R,x[1]=0,A[1]):
if nops(x)=1 then
 RETURN([seq(expand(coeff(R1,x[1],i)),i=0..A[1]-1)]):
else
[seq(INItay(coeff(R1,x[1],i),[op(2..n,x)],[op(2..n,A)]),i=0..A[1]-1)]:
fi:

end:



#Scheme(R,n,x,m,M): Let n be a positive integer (the second input parameter) and let R be a function of x[1],..., x[n]
#that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of
#the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m[1],...,m[n] 
#and the respective shift operators M[1],...,M[n]. Try:
#Scheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,m,M);
Scheme:=proc(R,n,x,m,M) local  i,j,lu,GU,INI:

GU:=[]:

for i from 1 to n do
lu:=MAZ(1,R/(mul(x[j]^(m[j]+1),j=1..i-1)*mul(x[j]^(m[j]+1),j=i+1..n)*x[i]),1/x[i],[seq(x[j],j=1..n)],m[i],M[i],{seq(m[j],j=1..i-1), seq(m[j],j=i+1..n)})[1]:
lu:=expand(1-subs(m[i]=m[i]-degree(lu,M[i]),lu)/M[i]^degree(lu,M[i])):
lu:=add(factor(coeff(lu,M[i],-j))*M[i]^(-j),j=1..-ldegree(lu,M[i])):
GU:=[op(GU),lu]:
od:

INI:=INItay(R,[seq(x[i],i=1..n)],[seq(-ldegree(GU[i],M[i]),i=1..n)]):
[GU,INI]:

end:



EvalScheme1:=proc(S,m,M,m0) local ope,INI,L,a,kha,j:
ope:=S[1][1]:
INI:=S[2]:
L:=INI:

if m0<=nops(INI)-1 then
 RETURN(INI[m0+1]):
fi:

for a from nops(INI) to m0 do

if member(0,
{seq(subs(m=a,denom(coeff(ope,M,-j))),j=1..-ldegree(ope,M))}
) then
 RETURN(FAIL):
else

kha:=normal(add(subs(m=a,coeff(ope,M,-j))*L[-j],j=1..-ldegree(ope,M))):
L:=[op(2..nops(L),L),kha]:
fi:
od:

L[-1]:

end:


#EvalScheme(S,m,M,m0): inputs a scheme S (let n:=nops(S)) in n dimension, phrased in terms of m and M (i.e. the discrete variables are m[1],...,m[n] and the corresponding shift operators are
#M[1],..., M[n], and a list of non-negative integers m0. Evaluats it at m0. Try
#EvalScheme([[2/M[1],3/M[2],4/M[3]],[[[1]]],m,M,[2,2,2]);
#EvalScheme([[m[1]/M[1]],[1]],m,M,[3]);

EvalScheme:=proc(S,m,M,m0) local n,ope1,GU,S1,A1,INI1,a,i1,m01,i,ka:

n:=nops(S[1]):
if nops(m0)<>n then
 RETURN(FAIL):
fi:

if n=1 then
RETURN(EvalScheme1(S,m[1],M[1],m0[1])):
fi:


ope1:=S[1][1]:
A1:=-ldegree(ope1,M[1]):




if m0[1]<A1 then
INI1:=S[2][m0[1]+1]:

S1:=[op(2..n,S[1])]:
S1:=subs(m[1]=m0[1],S1):
S1:=subs({seq(m[i1]=m[i1-1],i1=2..n),seq(M[i1]=M[i1-1],i1=2..n)},S1):
m01:=[op(2..n,m0)]:
ka:=EvalScheme([S1,INI1],m,M,m01):
if ka=FAIL then
RETURN(FAIL):
else
RETURN(EvalScheme([S1,INI1],m,M,m01)):
fi:
fi:

GU:=[]:

for a from 0 to A1-1 do

INI1:=S[2][a+1]:


S1:=[op(2..n,S[1])]:
S1:=subs(m[1]=a,S1):
S1:=subs({seq(m[i1]=m[i1-1],i1=2..n),seq(M[i1]=M[i1-1],i1=2..n)},S1):
m01:=[op(2..n,m0)]:

ka:=EvalScheme([S1,INI1],m,M,m01):

if ka=FAIL then
RETURN(FAIL):
else
GU:=[op(GU),ka]:
fi:
od:

ope1:=subs({seq(m[i]=m0[i],i=2..n)},ope1):

EvalScheme1([[ope1],GU],m[1],M[1],m0[1]):
end:



#CheckScheme(R,n,x,K1,K2): Checks SchemeS(R,n,x,m,M) for K2 different values in [0,K1]^n. Try:
#CheckScheme(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,10,10);
CheckScheme:=proc(R,n,x,K1,K2) local  S,m,M,ra,pt,i,i1:
S:=Scheme(R,n,x,m,M):
ra:=rand(0..K1):

for  i from 1 to K2 do

pt:=[seq(ra(),i1=1..n)]:

if EvalScheme(S,m,M,pt)<>CoefD(R,x,pt)  then
 print(`pt failed `):
 print( [EvalScheme(S,m,M,pt),CoefD(R,x,pt)]):
 RETURN(false):
fi:
od:

true:

end:




#SchemeV(R,n,x): A verbose version of Scheme(R,n,x,m,M) with illustration. try:
#SchemeV(1/(1-x[1]-x[2]-x[1]*x[2])^(1/3),2,x);
SchemeV:=proc(R,n,x) local  gu,m,M,F,i,t0,t1, ope,i1,i2,i3,INI,pt:

print(`A Pure Recurrence Scheme that enables Linear-Time and Constant-Space Calculation of the  Taylor coefficients of the function with`, n, `variables `, seq(x[i],i=1..n)):
print(``):
print(R):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, F(seq(m[i],i=1..n)), ` be the coefficient of`, mul(x[i]^m[i],i=1..n), ` in the Taylor expansion (around the origing) of the multivariable function (that may also be viewed as a formal power series) of`, seq(x[i],i=1..n) ):
print(``):
print(R):
print(``):

t0:=time():
gu:=Scheme(R,n,x,m,M):
t1:=time()-t0:

print(`The following PURE recurrence relations hold in each of the`, n, `discrete variables`, seq(m[i],i=1..n)):
print(``):

for i from 1 to n do
ope:=gu[1][i]:
print(``):
print(F(seq(m[i1],i1=1..n))= add(coeff(ope,M[i],-i1)*F(seq(m[i2],i2=1..i-1),m[i]-i1,seq(m[i2],i2=i+1..n)),i1=1..-ldegree(ope,M[i]) ) ):
print(``):
od:
print(``):
print(`By the way, it took`, t1, `seconds to geneate the scheme `):
print(``):
print(`Subject to the initial conditions`):
print(``):

INI:=gu[2]:

if n=1 then
 print(seq(F(i1)=INI[i1+1],i1=0..nops(INI)-1)):
elif n=2 then
 print(seq(seq(F(i1,i2)=INI[i1+1][i2+1],i2=0..nops(INI[i1+1])-1),i1=0..nops(INI)-1)):
elif n=3 then
 print(seq(seq(seq(F(i1,i2,i3)=INI[i1+1][i2+1][i3+1],i3=0..nops(INI[i1+1][i2+1])-1),i2=0..nops(INI[i1+1])-1),i1=0..nops(INI)-1)):
else
print(INI):
fi:

t0:=time():
print(`Using this scheme we can compute any value very fast. For example`):
pt:=[seq(1000*i1,i1=1..n)]:
print(``):
print(F(op(pt))=EvalScheme(gu,m,M,pt)):
print(``):
print(`This took`, time()-t0, `seconds. `):
print(``):

end:






#SchemeD(R,n,x,m,M): Let n be a positive integer (the second input parameter) and let R be a function of x[1],..., x[n]
#that has a "nice" Taylor expansion around the origin in n-space, uses the Apagodu-Zeilberger multivariable extension of
#the Almkvist-Zeilberger algorithm to output a PURE scheme phrased in terms of the discrete variables m
#and the  shift operator M for the coefficient of (x[1]*...*x[n])^m
#SchemeD(1/(1-x[1]-x[2]-x[3]-x[1]*x[2]*x[3])^(1/2),3,x,m,M);
SchemeD:=proc(R,n,x,m,M) local  i,j,lu,INI:



lu:=MAZ(1,R/mul(x[j],j=1..n),1/mul(x[j],j=1..n),[seq(x[j],j=1..n)],m,M,{})[1]:
lu:=expand(1-subs(m=m-degree(lu,M),lu)/M^degree(lu,M)):
lu:=add(factor(coeff(lu,M,-j))*M^(-j),j=1..-ldegree(lu,M)):
INI:=[seq(CoefD(R,x,[i$n]),i=0..-ldegree(lu,M)-1)]:
lu:=subs({M=M[1],m=m[1]},lu):
[[lu],INI]:

end:



#Walk2Dpaper(St,K): inputs a set of positive steps,St and outputs an article about the bi-sequence: Number of walks from the origin to a point [m1,m2].
#It gives the pure recurrence scheme  and computes as an illustrative example
#the number of ways of walking to [2*K,K] and [K,K]. It also gives the recurrence for the diagonal (i.e. the one-variable  and recurrence) for getting from the origin to [0,0] to [n,n]. For example, for
#the Royal walks : try:
#Walk2Dpaper({[1,0],[0,1],[1,1]},1000);
Walk2Dpaper:=proc(St,K) local i,s,R,x,Sc,a,b,A,B,m,M,OPES,INI,ope1,ope2,t0,lu,ope,i1,i2,Theorem:

t0:=time():
R:=1/(1-add(x[1]^s[1]*x[2]^s[2],s in St)):

Sc:=Scheme(R,2,x,m,M):
OPES:=Sc[1]:
INI:=Sc[2]:
OPES:=subs({m[1]=a,m[2]=b,M[1]=A,M[2]=B},OPES):

ope1:=OPES[1]:
ope2:=OPES[2]:

print(``):
print(` Theorem:  ` , `Let F(a,b) be the number of lattice walks, in the 2D Manhattan lattice,  from the origin to the point [a,b] using the atomic steps`):
print(``):
print(St):
print(``):
print(`Then F(a,b) satisfies the following two pure recurrences in the Eastbound and Northbound directions respectively`):
print(``):
print(F(a,b)=add(coeff(ope1,A,-i)*F(a-i,b),i=1..-ldegree(ope1,A))):
print(``):
print(F(a,b)=add(coeff(ope2,B,-i)*F(a,b-i),i=1..-ldegree(ope2,B))):
print(``):

print(`and in Maple notation`):
print(``):
lprint(F(a,b)=add(coeff(ope1,A,-i)*F(a-i,b),i=1..-ldegree(ope1,A))):
print(``):
lprint(F(a,b)=add(coeff(ope2,B,-i)*F(a,b-i),i=1..-ldegree(ope2,B))):
print(``):
print(`Subject to the initial conditions `):
print(``):
lprint(seq(seq(F(i1,i2)=INI[i1+1][i2+1],i1=0..nops(INI)-1),i2=0..nops(INI[1])-1)):
print(``):
print(`By the way, it took`, time()-t0, `seconds to generate the scheme `):




print(``):
t0:=time():
lu:=EvalScheme(Sc,m,M,[K,2*K]):
if lu<>FAIL and EvalScheme(Sc,m,M,[29,37])=Walk2Ddirect(St,[29,37]) then

print(`To illustrate how efficient it is, let's compute`, F(K,2*K), `and see how long it takes` ):
print(``):
print(F(K,2*K)=lu):
print(``):
print(`This took`, time()-t0, `seconds. Wow!`):
else
print(`There is a singularity, the scheme as it stands is useless`):
fi:

print(``):

t0:=time():
Sc:=SchemeD(R,2,x,m,M):
ope:=Sc[1][1]:
INI:=Sc[2]:

ope:=subs({m[1]=a,M[1]=A},ope):

print(``):
print(`Regarding the diagonal term, we have the following recurrence`):
print(``):
print(F(a,a)=add(coeff(ope,A,-i)*F(a-i,a-i),i=1..-ldegree(ope,A))):
print(``):

print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(F(i,i)=INI[i+1],i=0..nops(INI)-1)):

print(``):
print(`By the way, it took`, time()-t0, `seconds to generate this recurrence `):
print(``):
t0:=time():
print(`To illustrate how efficient it is, let's compute`, F(2*K,2*K), `and see how long it takes` ):
t0:=time():
lu:=EvalScheme(Sc,m,M,[2*K]):
print(``):
print(F(2*K,2*K)=lu):
print(``):
print(`This took`, time()-t0, `seconds. Wow!`):

if EvalScheme(Sc,m,M,[37])<>Walk2Ddirect(St,[37,37]) then
print(`Something bad happened, probably a singularity`):
 RETURN(FAIL):
fi:

print(``):
print(`Finally for the sake of Neil Sloane here are the first 30 terms of the diagonal sequence, starting at the walks to (1,1)`):
print(``):
print(seq(EvalScheme(Sc,m,M,[i]),i=1..30)):
print(``):
print(`-------------------------------------------------------------------`):
print(``):

end:



#Walk2DTh(St,K): inputs a set of positive steps,St and outputs an article about the bi-sequence: Number of walks from the origin to a point [m1,m2].
#It gives the pure recurrence scheme  and computes as an illustrative example
#the number of ways of walking to [2*K,K] and [K,K]. It also gives the recurrence for the diagonal (i.e. the one-variable  and recurrence) for getting from the origin to [0,0] to [n,n]. For example, for
#the Royal walks : try:
#Walk2DTh({[1,0],[0,1],[1,1]},1000);

Walk2DTh:=proc(St,K,nu) local i,s,R,x,Sc,a,b,A,B,m,M,OPES,INI,ope1,ope2,t0,lu,ope,i1,i2,Theorem:

t0:=time():
R:=1/(1-add(x[1]^s[1]*x[2]^s[2],s in St)):

Sc:=Scheme(R,2,x,m,M):
OPES:=Sc[1]:
INI:=Sc[2]:
OPES:=subs({m[1]=a,m[2]=b,M[1]=A,M[2]=B},OPES):

ope1:=OPES[1]:
ope2:=OPES[2]:

print(``):
print(cat( Theorem ,` ` , nu),  `: Let F(a,b) be the number of lattice walks, in the 2D Manhattan lattice,  from the origin to the point [a,b] using the atomic steps`):
print(``):
print(St):
print(``):
print(`Then F(a,b) satisfies the following two pure recurrences in the Eastbound and Northbound directions respectively`):
print(``):
print(F(a,b)=add(coeff(ope1,A,-i)*F(a-i,b),i=1..-ldegree(ope1,A))):
print(``):
print(F(a,b)=add(coeff(ope2,B,-i)*F(a,b-i),i=1..-ldegree(ope2,B))):
print(``):

print(`and in Maple notation`):
print(``):
lprint(F(a,b)=add(coeff(ope1,A,-i)*F(a-i,b),i=1..-ldegree(ope1,A))):
print(``):
lprint(F(a,b)=add(coeff(ope2,B,-i)*F(a,b-i),i=1..-ldegree(ope2,B))):
print(``):
print(`Subject to the initial conditions `):
print(``):
lprint(seq(seq(F(i1,i2)=INI[i1+1][i2+1],i1=0..nops(INI)-1),i2=0..nops(INI[1])-1)):
print(``):
print(`By the way, it took`, time()-t0, `seconds to generate the scheme `):




print(``):
t0:=time():
lu:=EvalScheme(Sc,m,M,[K,2*K]):
if lu<>FAIL and EvalScheme(Sc,m,M,[29,37])=Walk2Ddirect(St,[29,37]) then

print(`To illustrate how efficient it is, let's compute`, F(K,2*K), `and see how long it takes` ):
print(``):
print(F(K,2*K)=lu):
print(``):
print(`This took`, time()-t0, `seconds. Wow!`):
else
print(`There is a singularity, the scheme as it stands is useless`):
fi:

print(``):

t0:=time():
Sc:=SchemeD(R,2,x,m,M):
ope:=Sc[1][1]:
INI:=Sc[2]:

ope:=subs({m[1]=a,M[1]=A},ope):

print(``):
print(`Regarding the diagonal term, we have the following recurrence`):
print(``):
print(F(a,a)=add(coeff(ope,A,-i)*F(a-i,a-i),i=1..-ldegree(ope,A))):
print(``):

print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(F(i,i)=INI[i+1],i=0..nops(INI)-1)):

print(``):
print(`By the way, it took`, time()-t0, `seconds to generate this recurrence `):
print(``):
t0:=time():
print(`To illustrate how efficient it is, let's compute`, F(2*K,2*K), `and see how long it takes` ):
t0:=time():
lu:=EvalScheme(Sc,m,M,[2*K]):
print(``):
print(F(2*K,2*K)=lu):
print(``):
print(`This took`, time()-t0, `seconds. Wow!`):

if EvalScheme(Sc,m,M,[37])<>Walk2Ddirect(St,[37,37]) then
print(`Something bad happened, probably a singularity`):
 RETURN(FAIL):
fi:

print(``):
print(`Finally for the sake of Neil Sloane here are the first 30 terms of the diagonal sequence, starting at the walks to (1,1)`):
print(``):
print(seq(EvalScheme(Sc,m,M,[i]),i=1..30)):
print(``):
print(`-------------------------------------------------------------------`):
print(``):

end:


#BookWalks(N1,n1,K):  An article with lots of theorerms about lattice walks with sets of atomic stape of size from 2 to n1, and coordinates up to N1. 
#K is as in Walk2Dpaper(St,K) (q.v.) . Try:
#BookWalks(2,3,1000):

BookWalks:=proc(N1,n1,K) local W,i1,i2,n11,W1,St,nu,t0,St1:

print(`Efficient Counting of Lattice Walks For Many Choices of Sets of Allowed Steps in the 2-Dimensional Mahattan Lattice `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
t0:=time():

W:={seq(seq([i1,i2],i2=0..N1-i1),i1=0..N1)} minus {[0,0],[0,1],[1,0]}:

nu:=0:

for n11 from 0 to n1-2 do
W1:=choose(W,n11):

for St in W1 do

St1:=St union {[0,1],[1,0]}:

 nu:=nu+1:
print(``):
print(`----------------------------------------------------`):
print(``):
Walk2DTh(St1,K,nu):
print(``):

od:
od:

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

print(`This ends this book that took`, time()-t0, `seconds to produce `):

end:




#####NEEDS MORE WORK
#GJpaper(M,K): inputs a set of words M, in {1,2} and outputs a paper about the number of words in 1^a2^b avoiding the members of M as consectuive wubsord.. Try:
#GJpaper({[1,2,1],[2,1,2]},1000);

GJpaper:=proc(MIS,K) local i,R,x,Sc,a,A,m,M,t0,lu,ope,t,INI:

print(`Under construction still needs work`):
RETURN(FAIL):
t0:=time():
R:=normal(subs(t=0,GJgf({1,2},MIS,x,t))):

Sc:=SchemeD(R,2,x,m,M):

t0:=time():
Sc:=SchemeD(R,2,x,m,M):
ope:=Sc[1][1]:
INI:=Sc[2]:

ope:=subs({m[1]=a,M[1]=A},ope):

print(``):
print(`Theorem : Let F(a) be the number of words with a 1s and a 2s  that avoid, as consecutive subsords the following strings`):
print(``):
print(MIS):
print(``):
print(`Then F(a) satisfies the following recurrence `):
print(``):

print(``):
print(F(a)=add(coeff(ope,A,-i)*F(a-i),i=1..-ldegree(ope,A))):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(seq(F(i)=INI[i+1],i=0..nops(INI)-1)):
print(``):
print(`By the way, it took`, time()-t0, `seconds to generate this recurrence `):
print(``):
t0:=time():
print(`To illustrate how efficient it is, let's compute`, F(K), `and see how long it takes` ):
t0:=time():
lu:=EvalScheme(Sc,m,M,[K]):
print(``):
print(F(K)=lu):
print(``):
print(`This took`, time()-t0, `seconds. Wow!`):
print(``):
print(`Finally for the sake of Neil Sloane here are the first 30 terms of the diagonal sequence, starting at the walks to (1,1)`):
print(``):
print(seq(EvalScheme(Sc,m,M,[i]),i=1..30)):
print(``):
print(`-------------------------------------------------------------------`):
print(``):

end:



#####END NEEDS MORE WORK