with(combinat):
#######################################################################
## Save this file as  SMAZ.txt . Stay in the                          #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read `SMAZ.txt` <Enter>,                           #
## Then follow the instructions given there                           #.
##                                                                    #
## Written by Mohamud Mohammed and Doron Zeilberger,                  #
#Rutgers University ,                                                 #
##  [mohamudm, zeilberg] at math dot rutgers dot edu                  #
#######################################################################

with(combinat):
print(`First Written: Dec. 20, 2004: tested for Maple 8 and 9 `):
print(`Version of March 20 2011: (thanks for Miklos Bona) `):
print():
print(`This is SMAZ.txt, one of the Maple packages`):
print(`accompanying the article `):
#print(`Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and`):
#print(`the Sharpening of Wilf-Zeilberger Theory`):
#print(`by Moa Apagodu (formerly Mohammed Mohamud)  and Doron Zeilberger `):
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(`For general help, and a list of the available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
 
 

ezra:=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: