#DZtools.txt: For some of the projects in Math 640 (Rutgers Univ., Spring 2026, taught by Dr. Z.)

print(`DZtools.txt: A Maple package useful for Dr. Z.'s Final Projects for Math 640 (Experimental Mathematics, Rutgers Univeristy, Spring 2026`):
print(`For a list of the main procedures type: Help(); `):
print(`For a list of the other  procedures type: Help1(); `):
print(`With help with any procedures type: Help(ProcedureName);`):

Help:=proc():

if nargs=0 then
print(` The main procedures are: Findrec, RecDiaWalks, RecDiaWalksE, RecWalksBack, SeqFromRec `):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`It returns the guessed recurrence operator, followed by the initial conditions`):
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]=NuWalks then
print(`NuWalks(a,Steps): The number of walks from the origin to the point a in the d-dimens. lattice (where d:=nops(a)), using the set Steps all`):
print(`with non-neg. components, and at least one positive one, and of length d. Try:`):
print(`NuWalks([4,4,4],{[1,0,0],[0,1,0],[0,0,1]});`):

elif nargs=1 and args[1]=NuWalksDiaSeq then
print(`NuWalksDiaSeq(Steps,L): The first L terms of the sequence enumerating walks from the origin, using the set of steps, Steps, to the diagonal [n,..,n]. Try:`):
print(`NuWalksDiaSeq({[0,1],[1,0],[1,1]},30);`):

elif nargs=1 and args[1]=RecDiaWalks then
print(`RecDiaWalks(Steps,n,N): `):
print(`Given a set of symmetric steps (consisting of positive steps that are symmetric and of course the same dimension (size), let's call it d, uses`):
print(`the Mutlivariable Almkvist-Zeilberger algorithm to find a linear recurrence operator with polynomial coefficients in n, using  the positive`):
print(`shift operator N for the number of ways of walking from the origin to [n$d]. Try:`):
print(`RecDiaWalks({[1,0,0],[0,1,0],[0,0,1]}, n,N);`):

elif nargs=1 and args[1]=RecDiaWalksE then
print(`RecDiaWalksE(Steps,n,N,MaxC): `):
print(`Given a set of  steps (consisting of positive steps  uses`):
print(`non-rigorous (but very reliable) guessing to find an operator annihilating the sequence enumerating walks to [n,,,,.n] using steps`):
print(`If the set of steps, Steps, is symmetric, same output as RecDiaWalks(Steps,n,N); Try`):
print(`RecDiaWalksE({[1,0,0],[0,1,0],[0,0,1]}, n,N,4);`):

elif nargs=1 and args[1]=RecWalksBack then
print(`RecWalksBack(Steps,n,N): Given a set of symmetric steps (consisting of positive steps that are symmetric and of course the same dimension (size), let's call it d, uses`):
print(`the Mutlivariable Almkvist-Zeilberger algorithm to find a linear recurrence operator with polynomial coefficients in n, using  the set of steps and their reversals`):
print(`for the sequence enumerating the number of walks from the origin back to the origin using n steps. Try:`):
print(`RecWalksBack({[1,0,0],[0,1,0],[0,0,1]}, n,N);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(Ope,n,N,K): Given the pair Ope=[ope,Ini]`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec([N-n-1,[1]],n,N,10;`):

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

end:


Help1:=proc():

if nargs=0 then
print(` The supporting procedures are: NuWalks, NuWalksDiaSeq `):
else
Help(args):
fi:
end:

with(combinat):



#RecWalksBack(Steps,n,N): Given a set of symmetric steps (consisting of positive steps that are symmetric and of course the same dimension (size), let's call it d, uses
#the Mutlivariable Almkvist-Zeilberger algorithm to find a linear recurrence operator with polynomial coefficients in n, using  the set of steps and their reversals
#for the sequence enumerating the number of walks from the origin back to the origin using n steps. Try:
#RecWalksBack({[1,0,0],[0,1,0],[0,0,1]}, n,N);
RecWalksBack:=proc(Steps,n,N) local s,d,G,g,i1,x,f,y,z:
if not (type(Steps,set) and nops(Steps)>0 and nops({seq(nops(s), s in Steps)})=1 and min(seq(op(s), s in Steps))<=0) then
 RETURN(FAIL):
fi:
d:=nops(Steps[1]):


G:=permute(d):
for s in Steps do
 for g in G do
  if not member([seq(s[g[i1]],i1=1..d)],Steps) then
  print(`Steps are not symmetric`):
  RETURN(FAIL):
 fi: 
od:
od:

if not member([0$d],Steps) then

f:=add(mul(x[i1]^s[i1],i1=1..d), s in Steps)+add(mul(1/x[i1]^s[i1],i1=1..d), s in Steps):
else
f:=1+add(mul(x[i1]^s[i1],i1=1..d), s in Steps)+add(mul(1/x[i1]^s[i1],i1=1..d), s in Steps):
fi:


MAZ(1,1/mul(x[i1],i1=1..d),f,[seq(x[i1],i1=1..d)],n,N,{},[seq(y[i1],i1=1..d-1)],z)[1]:

end:

#RecDiaWalks(Steps,n,N): Given a set of symmetric steps (consisting of positive steps that are symmetric and of course the same dimension (size), let's call it d, uses
#the Mutlivariable Almkvist-Zeilberger algorithm to find a linear recurrence operator with polynomial coefficients in n, using  the positive
#shift operator N for the number of ways of walking from the origin to [n$d]. Try:
#RecDiaWalks({[1,0,0],[0,1,0],[0,0,1]}, n,N);
RecDiaWalks:=proc(Steps,n,N) local s,d,G,g,i1,x,f,y,z,ope:
if not (type(Steps,set) and nops(Steps)>0 and nops({seq(nops(s), s in Steps)})=1 and min(seq(op(s), s in Steps))<=0) then
 RETURN(FAIL):
fi:
d:=nops(Steps[1]):
if member([0$d],Steps) then
  RETURN(FAIL):
fi:

G:=permute(d):
for s in Steps do
 for g in G do
  if not member([seq(s[g[i1]],i1=1..d)],Steps) then
  print(`Steps are not symmetric`):
  RETURN(FAIL):
 fi: 
od:
od:

f:=1/(1-add(mul(x[i1]^s[i1],i1=1..d), s in Steps)):

f:=f/mul(x[i1],i1=1..d):

ope:=MAZ(1,f,1/mul(x[i1],i1=1..d),[seq(x[i1],i1=1..d)],n,N,{},[seq(y[i1],i1=1..d-1)],z)[1]:

[ope,NuWalksDiaSeq(Steps,degree(ope,N))]:
end:

 

#FROM SMAZ.txt
ezraSMAZ:=proc()
if nargs=0 then
print(`SMAZ.txt :`):
print(`A Maple package for  finding recurrences for multi-integrals of `):
print(` Hypergeometric/Hyperexponential with symmetric Integrands.`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`The main procedure is:  DIAGpaper, MAZ,MAZm  `):

elif nargs=1 and args[1]=DIAGpaper then
print(`DIAGpaper(TOP,BOT,a,x,n,pars,F,A): Inputs polynomials in the list of variables`):
print(`x, TOP, and BOT, and a numeric or symbolic number a, and a symbol m, and a set of parameters pars`):
print(`(if you are lazy, you can make it {}) and letters F,A,R, (for formulating the proof)`):
print(`outputs a mathematical article`):
print(`authored by Shalosh B. Ekhad that finds and proves a linear recurrence equation for the`):
print(`diagonal sequence (called A(n)), i.e. the coefficient of x[1]^n*x[2]^n* in the formal power series`):
print(`TOP/BOT^a. BOT must have a non-zero constant terms (for the formal power series to be well-defined)`):
print(`For example, for Miklos Bona's problem, try:`):
print(`DIAGpaper(1,(1-x)^2+(1-y)^2-1,1/2,[x,y],n,{},F,A);`):
 
elif nargs=1 and args[1]=MAZ then
print(`MAZ(POL,H,rat,x,n,N,pars,y,z): Given a polynomial, POL, an hyper-exponential function, H, and`):
print(`a rational function, rat`):
print(`of the continuous variables in the list x, all of which are SYMMETRIC w.r.t. the variables in [x]`):
print(`and a discrete variable`):
print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, pars`):
print(`ouputs a recurrence operator, let's call it  ope(n,n), and a multi-certificate, let's call it R(z;y)`):
print(`such that the function F(n;x):=POL*H*rat^n satisfies`):
print(` ope(n,N)F=sum(D_{x[i]} [(R(x[i];x[1], ..., x[i-1],x[i+1]..)H*rat^n)],i=1..nops(x)`):
print(`In particular, try:`):
print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{},[y1],z1);`):
print(`MAZ(1,exp(-x^2/2-y^2/2-z^2/2),((x-y)*(x-z)*(y-z))^2,[x,y,z],n,N,{},[y1,y2],z1);`):

print(`MAZ(x+y+z,1/(1-x-y-z+a*x*y*z)/x/y/z,1/x/y/z,[x,y,z],n,N,{a},[y1,y2],z1);`):


elif nargs=1 and args[1]=MAZm then
print(`MAZm(POL,H,rat,x,n,N,pars,y,z): Like MAZ but with the lsit y`):
print(`replaced by the letter m, and the second output, R,`):
print(` is given in terms of`):
print(`monomial symmetric functions using m`):

print(`For example, try:`):
print(`MAZm(1,exp(-x^2/2-y^2/2-z^2/2),((x-y)*(x-z)*(y-z))^2,[x,y,z],n,N,{},m,Z);`):


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

fi:



end:




####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,y,z): 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,y,z)
local deg,gu,i,i1:
deg:=FindDeg(POL,H,rat,x[1],{seq(x[i1],i1=1..nops(x))},n,L):

for i from 0 to deg do
 gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,i,y,z):

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


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

for i from 1 to nops(x)-1 do
if expand(subs({x[i]=x[i+1],x[i+1]=x[i]},POL)-POL)<>0
or normal(subs({x[i]=x[i+1],x[i+1]=x[i]},H)/H)<>1
or normal(subs({x[i]=x[i+1],x[i+1]=x[i]},rat)/rat)<>1
  then 
   ERROR(`Integrand not symmetric in the variables in`,x):
fi:
od:

for L from 0 do
#print(`trying L=`,L):
gu:=MAZ1(POL,H,rat,x,n,N,L,pars,y,z):

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


CheckMAZ:=proc(POL,H,rat,x,n,N,y,z,ope,R)
local F,luL,luR,i,Rlist,i1:
F:=H*rat^n:
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))):

Rlist:=[]:

for i from 1 to nops(x) do
Rlist:=
[op(Rlist),subs({z=x[i],seq(y[i1]=x[i1],i1=1..i-1),seq(y[i1]=x[i1+1],
i1=i..nops(x)-1)},R)]:
od:

luR:=normal(add( normal(diff(Rlist[i]*F,x[i])/F ),
              i=1..nops(Rlist))):

evalb(normal(luR-luL)=0):
end:





#MAZ1deg(POL,H,rat,x,n,N,L,pars,deg,y,z): 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,y,z)
local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,i1,eqN,var1N,opeN,bilti,meka,
R,ope1,R1:

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:

X:=[]:

R:=GenASP(deg,z,y,a):
var:=R[2]:
R:=R[1]:

for i from 1 to nops(x) do
 X:=[op(X),subs({z=x[i],seq(y[i1]=x[i1],i1=1..i-1),
seq(y[i1]=x[i1+1],i1=i..nops(y))}, R)]:
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):
else
# print(`There is hope for a recurrence of order`,L):
fi:

var1:=solve(eq,var):

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

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
ope1:=coeff(ope,bilti[i],1):
if ope1<>0 then
  R1:=coeff(R,bilti[i],1):
  meka:=coeff(ope1,N,degree(ope1,N)):
  ope1:=ope1/meka: R1:=normal(R1/meka):
  RETURN(yafe(ope1,N),R1):
 fi:
od:

end:


####About Symmetric Polynomials
ezra1:=proc()
if args=NULL then
print(`The supportingprocedures are:  CheckMAZ`):
print(`GenASP, GenSP,mlam, Pars, Pars1 , `):
print():
fi:


if nargs=1 and args[1]=GenASP then 
print(`GenASP(n,z,y,a): An almost-symmetric generic polynomial with distinguished variable`):
print(`z and symmetric variables y of degree n, indexed by a, followed by the`):
print(`set of coefficients`):
fi:

if nargs=1 and args[1]=GenSP then 
print(`GenSP(n,x,a): a generic symmetric polynomial of degree n in the varialbes x`):
print(`indexed by a, followed by the set of coefficients`):
fi:

if nargs=1 and args[1]=mlam then 
print(`mlam(lam,x): given a partition lam and a list of variables of the`):
print(`same length, finds the monomial symmetric function m_lam(x)`):
fi:

if nargs=1 and args[1]=Pars then
print(`Pars(n,k): all the partitions of n into k parts (including 0)]`):
fi:

if nargs=1 and args[1]=Pars1 then
print(`Pars1(n,k,a): partitions of n into exactly k parts `):
print(`whose largest part is a`):
fi:


end:


#Pars1(n,k,a): partitions of n into exactly k parts 
#whose largest part is a
Pars1:=proc(n,k,a)
local gu,i,a1,gu1,n1:

if a=0 then
 if n=0 then
   RETURN({[0$k]}):
 else
   RETURN({}):
 fi:
fi:

if k=1 then
 if n=a then
  RETURN({[a]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:

n1:=n-a:

for a1 from 0 to a do
 gu1:=Pars1(n1,k-1,a1):
 gu:=gu union  {seq([a,op(gu1[i])],i=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:

#Pars(n,k): all the partitions of n into k parts (including 0)]
Pars:=proc(n,k)
local a1:
{seq(op(Pars1(n,k,a1)),a1=0..n)}:
end:


#mlam(lam,x): given a partition lam and a list of variables of the
#same length, finds the monomial symmetric function m_lam(x)
mlam:=proc(lam,x) local i,gu,j:

if nops(x)<>nops(lam) then
  ERROR(`Both inputs must be lists of the same length`):
fi:

gu:=permute(lam):
add(mul(x[j]^gu[i][j],j=1..nops(x)),i=1..nops(gu)):
end:


#GenSP(n,x): a generic symmetric polynomial of degree n in the varialbes x
#followed by the set of coefficients
GenSP:=proc(n,x,a)
local gu,k,n1,i,kv,lu:
k:=nops(x):

lu:={seq(op(Pars(n1,k)),n1=0..n)}:
kv:={}:

gu:=0:

for i from 1 to nops(lu) do
gu:=gu+a[i]*mlam(lu[i],x):
kv:=kv union {a[i]}:
od:
gu,kv:
end:


#GenASP(n,z,y,a): An almost-symmetric generic polynomial with distinguished variable
#z and symmetric variables y of degree n, indexed by a, followed by the
#set of coefficients
GenASP:=proc(n,z,y,a)
local gu,n1,kv,mu,i:
if y=[] then
 RETURN(add(a[i]*z^i,i=0..n),{seq(a[i],i=0..n)}):
fi:
kv:={}:
gu:=0:
for n1 from 0 to n do
mu:=GenSP(n-n1,y,a[n1]):
gu:=gu+z^n1*mu[1]:
 kv:=kv union mu[2]:
od:
gu,kv:
end:

########End Symmetric Polynomials


#SPtom(P,y,m): Given a symmetric polynomial, outputs
#it in terms of the monomial symmetric functions

SPtom:=proc(P,x,m)
local i,gu,mu,P1,NewP,PLO:

if x=[] then
  RETURN(P):
fi:

for i from 1 to nops(x)-1 do
if expand(subs({x[i]=x[i+1],x[i+1]=x[i]},P)-P)<>0
  then 
   ERROR(`P is not symmetric in the variables in`,x):
fi:
od:

if not type(P,`+`) then
     gu:=[seq(degree(P,x[i]),i=1..nops(x))]:
     if nops({op(gu)})<>1 then
           ERROR(`Not symmetric`):
     fi:
     mu:=normal(P/mul(x[i]^gu[i],i=1..nops(x))):
      RETURN(mu*m[op(gu)]):
fi:

NewP:=0:
PLO:=P:

while PLO<>0 do
 P1:=op(1,PLO):
     gu:=[seq(degree(P1,x[i]),i=1..nops(x))]:
     mu:=PLO:
        for i from 1 to nops(x) do
           mu:=coeff(mu,x[i],gu[i]):
        od:  
     gu:=sort(gu,`>`):
     NewP:=NewP+factor(mu)*m[op(gu)]:
     PLO:=expand(PLO-mu*mlam(gu,x)):
od:
NewP:
     
end:



# MAZm(POL,H,rat,x,n,N,pars,m,z):Like MAZ but the certificate in m notation
#z is still the main variable
MAZm:=proc(POL,H,rat,x,n,N,pars,m,z) local gu,i,ope,R,Y,y,Rtop,Rbot:

Y:=[seq(y[i],i=1..nops(x)-1)]:

gu:=MAZ(POL,H,rat,x,n,N,pars,Y,z):



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

ope:=gu[1]:

R:=gu[2]:
R:=normal(R):
Rtop:=numer(R):
Rbot:=denom(R):

Rtop:=add(SPtom(coeff(Rtop,z,i),Y,m)*z^i,i=0..degree(Rtop,z)):
ope,Rtop/Rbot:

end:

#DIAGpaper(TOP,BOT,a,x,n,pars,F,A): Inputs polynomials in the list of variables
#x, TOP, and BOT, and a numeric or symbolic number a, and a symbol m, and a set of parameters pars
#(if you are lazy, you can make it {}) and letters F,A,R, (for formulating the proof)
#outputs a mathematical article
#authored by Shalosh B. Ekhad that finds and proves a linear recurrence equation for the
#diagonal sequence (called A(n)), i.e. the coefficient of x[1]^n*x[2]^n* in the formal power series
#TOP/BOT^a. BOT must have a non-zero constant terms (for the formal power series to be well-defined)
#For example, for Miklos Bona's problem, try:
#DIAGpaper(1,(1-x)^2+(1-y)^2-1,1/2,[x,y],n,{},F,A):

DIAGpaper:=proc(TOP,BOT,a,x,n,pars,F,A) local gu,i,j,ope,N,t0,y,z,R:

t0:=time():
if subs({seq(x[i]=0,i=1..nops(x))},BOT)=0 then
 print(`The bottom`, BOT, `has zero constant-terms , so the formal power series is undefined`):
 RETURN(FAIL):
fi:


gu:=MAZ(TOP,1/BOT^a/mul(x[i],i=1..nops(x)),1/mul(x[i],i=1..nops(x)),x,n,N,pars,[seq(y[j],j=1..nops(x)-1)],z):

print(``):

ope:=gu[1]:

if gu=FAIL then
 RETURN(FAIL):
fi:
print(`A Linear Recurrence Equation for The Diagonal Coefficients of The  Power Series Of`):
print(TOP/BOT^a):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` Let `, A(n), `be the coefficient of`, mul(x[i]^n,i=1..nops(x)), `in the Maclaurin Expansion of the Formal Power series`,  TOP/BOT^a ):
print(`then `, A(n), `satisfies the following linear recurrence equation with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Proof:`):
print(`Thanks to Cauchy's Integral Formula, A(n) equals to  `):
print(` a constant (independent of n) times `):
print(`the contour-integral with respect to the complex variables`, op(x)):
print(`around any poly-circle around the origin, of the function`):
print(F(n,op(x))=TOP/BOT^a/mul(x[i]^(n+1),i=1..nops(x))):
print(`Let's cleverly constuct rational function`):

print( R(seq(y[j],j=1..nops(x)-1),z)= gu[2]):
print(`that certifies it`):
print(``):
print(`(Check!)`):
print(`and the theorem follows upon contour-integrating with respect to`, op(x)):
print(` QED! `):

print(`This took`, time()-t0, `seconds . `):

end:




###From FindRec.txt
######################################################################
## FindRec.txt: Save this file as FindRec.txt To use it, stay in the #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read FindRec.txt <Enter>                          #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  zeilberg@math.rutgers.edu.                                       # 
######################################################################


print(`First Written: April 15, 2005: tested for Maple 8 and 9 `):
print(`Version of April 15, 2005:  `):
print():
print(`This is FindRec.txt A Maple package`):
print(`accompanying Doron Zeilberger's  article: `):
print(`The Holonomic Ansatz: I. Foundations`):
print(`It guesses recurrence equations for discrete variables of ONE variable`):
print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
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():
 
 

 
ezraFindRec:=proc()
if args=NULL then

print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


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


fi:
 
end:


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

ope:=subs(var1,ope):

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

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

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

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

end:





 
#SeqFromRecOld(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRecOld:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

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

gu:

end:

#SeqFromRec(Ope,n,N,K): Given the pair Ope=[ope,Ini]
#extends it to the first K values
SeqFromRec:=proc(Ope,n,N,K)
local ope,Ini,ope1,gu,L,n1,j1:
ope:=Ope[1]:
Ini:=Ope[2]:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

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

gu:

end:

#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
###End From FindRec.txt



#NuWalks(a,Steps): The number of walks from the origin to the point a in the d-dimens. lattice using the set Steps. Try:
#NuWalks([4,4,4],{[1,0,0],[0,1,0],[0,0,1]});
NuWalks:=proc(a,Steps) local d,s:
option remember:
d:=nops(a):
if not (type(a,list)  and {seq(nops(s), s in Steps)}={d} ) then
 RETURN(FAIL):
fi:

if a=[0$d] then
 RETURN(1):
fi:
if min(a)<0 then
 RETURN(0):
fi:
add(NuWalks(a-s,Steps), s in Steps):
end:


#NuWalksDiaSeq(Steps,L): The first L terms of the sequence enumerating walks from the origin, using the set of steps, Steps, to the diagonal [n,..,n]. Try:
#NuWalksDiaSeq({[0,1],[1,0],[1,1]},30);
NuWalksDiaSeq:=proc(Steps,L) local s,d,n:
if  nops({seq(nops(s), s in Steps)})<>1 and type(L, integer) and L>0  then
 RETURN(FAIL):
fi:

d:=nops(Steps[1]):
[seq(NuWalks([n$d],Steps),n=1..L)]:

end:




#RecDiaWalksE(Steps,n,N,MaxC): same output as RecDiaWalks(Setps,n,N) but non-rigorously by guessing, trying to find an operator of complexity MaxC. Try:
#RecDiaWalksE({[1,0,0],[0,1,0],[0,0,1]}, n,N,4);
RecDiaWalksE:=proc(Steps,n,N,MaxC) local L,ORDER:

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

Findrec(NuWalksDiaSeq(Steps,L),n,N,MaxC):

end: