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


print(`First Written: Dec. 7, 2004: tested for Maple 8 and 9 `):
print(`Version of Dec. 7, 2004:  `):
print():
print(`This is MultiAlmkvistZeilberger, 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 Mohamud Mohammed 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(`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 `):
 
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,{});`):


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

fi:



end:



ez:=proc():print(`MAZ1(POL,H,rat,x,n,N,L,pars), MAZ(POL,H,rat,x,n,N,pars)`):
print(`FindDeg(POL,H,rat,x,n,L)`):
print(`CheckMAZ1(POL,H,rat,x,n,N,L,pars)`):
end:


_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
print(`trying L=`,L):
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:



#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):
else
  print(`There is hope for a recurrence of order`,L):
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:


if nops(bilti)>1 then
 print(`There is some slack`):
fi:

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: