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

#Created: June 4, 2002
#This version: June 4, 2002
#WZtriple: A Maple Program for finding a linear recurrence equation with
#polynomial coefficients, w.r.t. n, for expressions of the form
#Int_{z,w,x}F(n,z,w,x), where F is a hypergeometric term
#See Wilf and Zeilberger, Invent. math. 108(1992) 575-633
#A better version is available by Akalu Tefera
#Please report bugs to zeilberg@math.rutgers.edu




print(`Created: June 4, 2002`):
print(`This version: June 4, 2002`):
print(`WZtriple: A Maple Program for finding a linear recurrence`):
print(` equation with polynomial coefficients, w.r.t. n, `):
print(`for expressions of the form Int_{z,w,x}F(n,z,w,x), where F `):
print(`is a hypergeometric term`):
print(`It is based on the method in Wilf's and Zeilberger's 1992 article`):
print(` An algorithmic proof theory for hypergeometric`):
print(` (ordinary and "q") multisum/integral identities}, Invent. Math. `):
print(` v. 108, 575-633 (1992). Available from Wilf's and Zeilberger's`):
print(` websites`):
print(``):
print(` A better, and more general version was made by Akalu Tefera`):
lprint(``):
lprint(` IMPORTANT WARNING: a, R1,R2,R3,Y1R1,Y2R2,Y3R3 are global `):
print(` variables, do not use them! `):

lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.rutgers.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):
ezra:=proc()
if args=NULL then
 print(`Contains the following procedures: WZ2, WZ3, WZ3s `):
fi:

if nops([args])=1 and op(1,[args])=WZ2 then
print(`WZ2(F,n,x,y,resh,N):`):
print(`Given a  closed-form function F`):
print(`in the discrete variable n and the continuous variables x and y`):
print(`tries to find a recurrence operator ope(n,N)`):
print(`(where N is the shift in n (Nf(n)=f(n+1)),`):
print(` ope(N,n)F(n,x,y)=D_x(G1*F)+D_y(G2*F) `):
print(`It returns: ope,G1,G2 in case of success`):
print(` 0 otherwise `):
print(``):
print(`For example, type:`):
print(` WZ2((x+y)^(2*n)/x^(n+1)/y^(n+1),n,x,y,[n],N):`):

fi:


if nops([args])=1 and op(1,[args])=WZ3 then
print(`WZ3(F,n,x,y,z,resh,N):`):
print(`Given a  closed-form function F`):
print(`in the discrete variable n and the continuous variables x, y, and z`):
print(`tries to find a recurrence operator ope(n,N)`):
print(`(where N is the shift in n (Nf(n)=f(n+1)),`):
print(` ope(N,n)F(n,x,y,z)=D_x(G1*F)+D_y(G2*F)+D_z(G3*F) .`):
print(`It returns: ope,G1,G2,G3 in case of success`):
print(` 0 otherwise `):
print(``):
print(`For example, type:`):
print(` WZ3((x+y+z)^(3*n)/x^(n+1)/y^(n+1)/z^(n+1),n,x,y,z,[n],N):`):

fi:

if nops([args])=1 and op(1,[args])=WZ3s then
print(`WZ3s(P,F,n,x,y,z,ORDER,resh,N,S1,S2,S3):`):
print(`Given a polynomial P and closed-form function F`):
print(`in the discrete variable n and the continuous variables`):
print(`x,y,z, and a non-negative integer ORDER, and rational functions`):
print(`(in x,y,z), S1,S2,S3 `):
print(`tries to find a recurrence operator ope(n,N)`):
print(`(where N is the shift in n (Nf(n)=f(n+1)), of order (degree in N))`):
print(` ORDER, and certificates G1,G2,G3, such that `):
print(` ope(N,n)F(n,x,y,z)=D_x(G1*S1*F)+D_y(G2*S2*F)+D_z(G3*S3*F) `):
print(`It returns: ope,G1,G2,G3 in case of success`):
print(` 0 otherwise `):
print(``):
print(`For example, type:`):
print(`P:=1,F:=(x+y+z)^(3*n)/x^(n+1)/y^(n+1)/z^(n+1):ORDER:=1:`):
print(` resh:=[n]: , S1:=1/x/y/z: S2:=1/x/y/z: S3:=1/x/y/z:`):
print(` WZ3s(P,F,n,x,y,z,ORDER,resh,N,S1,S2,S3):`):

fi:

end:






#findden:Given a differential equation Left-coeffR1*R1-D_y1R1-
#coeffR2*R2-D_y2R2 finds the "generic" denominators
#that R1,R2 should have

#cv is  a program that inputs a functional equation,eq, in
#R1(y1,y2), Y1R1(y1,y2)(D_y1 R1),R2(y1,y2),Y2R2(y1,y2)(=D_y2 R2)
#and some proposed Change of dependent Variables(and hence called cv)
#OldR1:=NewR1*S1, OldR2:=NewR2*S2. The output is the functional equation
#in the NewR1,NewR2. Another part of the input is S1P (P for previous)
#S2P, where R1Ancient=S1P*OldR1,R2Ancient=S1P*OldR2, the second and third
#part of the output is the realtion between R1Ancient and R1new, etc,
#i.e. S1P*S1, S2P*S2

findden:=proc(eq,R1)
local Left,coeffR1,gu:
Left:=coeff(eq,Smol,1):
coeffR1:=coeff(eq,R1,1):
Left:=normal(Left):
gu:=denom(coeffR1)/denom(Left):
gu:=normal(gu):
1/denom(gu):
end:

pashet:=proc(p,N)
local i,gu,p1:
p1:=expand(p):
gu:=0:
for i from 0 to degree1(p,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu):
end:

nor:=proc(bitu)
 local i,gu:
gu:=0:

if type(bitu,`+`) then

 for i from 1 to nops(bitu) do 
  gu:=gu+normal(op(i,bitu)):
 od:

else
 gu:=normal(bitu):

 fi:
 gu:

end:



cv:=proc(S1P,S2P,S3P,S1,S2,S3,eq,y1,y2,y3)
local gu:
gu:=eq:
gu:=subs(R1=R1*S1,gu):
gu:=subs(R2=R2*S2,gu):
gu:=subs(R3=R3*S3,gu):
gu:=subs(Y1R1=S1*Y1R1+diff(S1,y1)*R1,gu):
gu:=subs(Y2R2=S2*Y2R2+diff(S2,y2)*R2,gu):
gu:=subs(Y3R3=S3*Y3R3+diff(S3,y3)*R3,gu):
nor(gu),simplify(S1*S1P),simplify(S2*S2P),simplify(S3*S3P):
end:




#This program fa, finds the functional equation for the rational functions
#R1(n,y1,y2),R2(n,y1,y2), satisfying, for a given CF 
#P(n,y1,y2) F(n,y1,y2) (P(n,y1,y2) being polynomial) the functional
#equation 
#Operator(N,n)F
#D_y1 (R1(n,y1,y2)F)+
#D_y2 (R2(n,y1,y2)F(n,y1,y2))

#rf converts (a)_y (Raising Factorial) to GAMMA
rf:=proc(a,y)

GAMMA(a+y)/GAMMA(a):

end:

fa:=proc(P,F,n,y1,y2,y3,ORDER,R1,R2,R3,Y1R1,Y2R2,Y3R3)
local i,gu:

gu:=P:

for i from 1 to ORDER do
 gu:=gu+a[i]*subs(n=n+i,P)*simplify(subs(n=n+i,F)/F):
od:

gu:=gu*Smol- Y1R1-R1*normal(diff(simplify(log(F)),y1))
       - Y2R2-R2*normal(diff(simplify(log(F)),y2))
-Y3R3-R3*normal(diff(simplify(log(F)),y3)):

end:










#A program that multiplies a polynomial pol by an equation
#in R1,R2,Y1R1,Y2R2


finddeg:=proc(eq,R1,Y1R1,x):
degree1(eq,x)-max(degree1(coeff(eq,R1),x),degree1(coeff(eq,Y1R1),x)):
end:


multpa:=proc(eq,pol)
local pu,i:


if type(eq,`+`) then

 pu:=0:

  for i from 1 to nops(eq) do
   pu:=pu+normal(op(i,eq)*pol):
  od:
  
 RETURN(pu):

else

 RETURN(normal(eq*pol)):

fi:


end:






#A program that given a differential equation in R1(y1,y2),
#Y1R1:=D_y1 R1(y1,y2), R2(y1,y2), Y2R2:=D_y2 R2(y1,y2), and
#prescribed degrees degR1y1,degR1y2,tdegR1,degR2y1,degR2y2,tdegR2
#and prescribed ORDER
#solves a given equation eq=0
#ptorpols(eq,y1,y2,degR1y1,degR1y2,tdegR1,degR2y1,degR2y2,tdegR2)
#eq is the equation,y1,y2, are the variables
#Procedure genpol generates a generic polynomial with degree dk in k 
#and degree dl in l, and total degree tdeg, the coefficients are
#a[i,j]

spec:=proc(a,dk,dl,dm,tdeg,davar)
local i,j,j1,gu:
gu:=davar:
for i from 0 to dk do
for j from 0 to min(dl,tdeg-i) do
for j1 from 0 to min(dm,tdeg-i-j) do
gu:=subs(a[i,j,j1]=0,gu):
od:
od:
od:
RETURN(gu):
end:




genpol:=proc(a,k,l,m,dk,dl,dm,tdeg)
local i,j,j1,gu,gu1:
gu1:={}:
gu:=0:
for i from 0 to dk do
for j from 0 to min(dl,tdeg-i) do
for j1 from 0 to min(dm,tdeg-i-j) do
gu:=gu+a[i,j,j1]*k^i*l^j*m^j1:
gu1:=gu1 union {a[i,j,j1]}:
od:
od:
od:
RETURN(gu,gu1):
end:


ptorpols:=proc(e,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3,ORDER,N)
local i,j,j1,eq,equ,var,pol1,pol2,pol3,Y1pol1,Y2pol2,
Y3pol3,d1,d2,d3,d4,glu,lu,ope:

ope:=1:

for i from 1 to ORDER do
 ope:=ope+a[i]*N^i:
 od:

var:={}:

for i from 1 to ORDER do
 var:=var union {a[i]}
od:

eq:=e:
pol1:=genpol(a,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1):
var:=var union pol1[2]:
pol1:=pol1[1]:

pol2:=genpol(b,y1,y2,y3,degR2y1,degR2y2,degR2y3,tdegR2):
var:= var union pol2[2]:
pol2:=pol2[1]:

pol3:=genpol(c,y1,y2,y3,degR3y1,degR3y2,degR3y3,tdegR3):
var:= var union pol3[2]:
pol3:=pol3[1]:

Y1pol1:=diff(pol1,y1):
Y2pol2:=diff(pol2,y2):
Y3pol3:=diff(pol3,y3):

eq:=subs(R1=pol1,eq):
eq:=subs(R2=pol2,eq):
eq:=subs(R3=pol3,eq):
eq:=subs(Y1R1=Y1pol1,eq):
eq:=subs(Y2R2=Y2pol2,eq):
eq:=subs(Y3R3=Y3pol3,eq):
eq:=expand(eq):

equ:={}:
d1:=degree1(eq,y1):
d2:=degree1(eq,y2):
d3:=degree1(eq,y3):
d4:=degree1(eq,{y1,y2,y3}):

for i from 0 to d1 do
 for j from 0 to min(d2,d4-i) do
  for j1 from 0 to min(d3,d4-i-j) do
  glu:=coeff(coeff(coeff(eq,y1,i),y2,j),y3,j1):
   if glu<>0 then
    equ:=equ union {glu=0}:
   fi:
  od:
 od:
od:


lu:=solve(equ,var):

if lu=NULL then
 RETURN(0):
fi:

pol1:=subs(lu,pol1):
pol2:=subs(lu,pol2):
pol3:=subs(lu,pol3):
ope:=subs(lu,ope):

pol1:=spec(a,degR1y1,degR1y2,degR1y3,tdegR1,pol1):
pol1:=spec(b,degR2y1,degR2y2,degR2y3,tdegR2,pol1):
pol1:=spec(c,degR3y1,degR3y2,degR3y3,tdegR3,pol1):

pol2:=spec(a,degR1y1,degR1y2,degR1y3,tdegR1,pol2):
pol2:=spec(b,degR2y1,degR2y2,degR2y3,tdegR2,pol2):
pol2:=spec(c,degR3y1,degR3y2,degR3y3,tdegR3,pol2):


pol3:=spec(a,degR1y1,degR1y2,degR1y3,tdegR1,pol3):
pol3:=spec(b,degR2y1,degR2y2,degR2y3,tdegR2,pol3):
pol3:=spec(c,degR3y1,degR3y2,degR3y3,tdegR3,pol3):

ope:=spec(a,degR1y1,degR1y2,degR1y3,tdegR1,ope):
ope:=spec(b,degR2y1,degR2y2,degR2y3,tdegR2,ope):
ope:=spec(c,degR3y1,degR3y2,degR3y3,tdegR3,ope):
pol1:=factor(pol1):
pol2:=factor(pol2):
pol3:=factor(pol3):
RETURN(ope,pol1,pol2,pol3);

end:




#Like ptorpolstik, but with random assignements of the parameters in
#the list resh

ptorpolstik:=proc(e,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3, ORDER,resh)

local i,j,j1,eq,equ,var,pol1,pol2,pol3,Y1pol1,Y2pol2,Y3pol3,
d1,d2,d3,d4,glu,lu,ra:

ra:=2:

eq:=e:
for i from 1 to nops(resh) do
 eq:=subs(op(i,resh)=ra,eq):
od:

var:={}:

for i from 1 to ORDER do
 var:=var union {a[i]}
od:


pol1:=genpol(a,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1):
var:=var union pol1[2]:
pol1:=pol1[1]:

pol2:=genpol(b,y1,y2,y3,degR2y1,degR2y2,degR2y3,tdegR2):
var:= var union pol2[2]:
pol2:=pol2[1]:

pol3:=genpol(c,y1,y2,y3,degR3y1,degR3y2,degR3y3,tdegR3):
var:= var union pol3[2]:
pol3:=pol3[1]:


Y1pol1:=diff(pol1,y1):
Y2pol2:=diff(pol2,y2):
Y3pol3:=diff(pol3,y3):
eq:=subs(R1=pol1,eq):
eq:=subs(R2=pol2,eq):
eq:=subs(R3=pol3,eq):
eq:=subs(Y1R1=Y1pol1,eq):
eq:=subs(Y2R2=Y2pol2,eq):
eq:=subs(Y3R3=Y3pol3,eq):
eq:=expand(eq):

equ:={}:
d1:=degree1(eq,y1):
d2:=degree1(eq,y2):
d3:=degree1(eq,y3):
d4:=degree1(eq,{y1,y2,y3}):

for i from 0 to d1 do
 for j from 0 to min(d2,d4-i) do
  for j1 from 0 to min(d3,d4-i-j) do
  glu:=coeff(coeff(coeff(eq,y1,i),y2,j),y3,j1):
   if glu<>0 then
    equ:=equ union {glu=0}:
   fi:
  od:
 od:
od:



lu:=solve(equ,var):


if lu<> NULL then
 RETURN(1):
else
 RETURN(0):
fi:


end:

WZ3s:=proc(P,F,n,y1,y2,y3,ORDER,resh,N,S1,S2,S3)
local eq,gu,degR1y1,degR1y2,degR1y3,tdegR1,degR2y1,degR2y2,degR2y3,tdegR2,
degR3y1,degR3y2,degR3y3,tdegR3,ope,R_1,R_2,R_3,mekh:
eq:=fa(P,F,n,y1,y2,y3,ORDER,R1,R2,R3,Y1R1,Y2R2,Y3R3):
#print(`eq is`,eq):
eq:=subs(Smol=1,eq):
gu:=cv(1,1,1,S1,S2,S3,eq,y1,y2,y3):
eq:=gu[1];
eq:=numer(normal(eq)):


degR1y1:=finddeg(eq,R1,Y1R1,y1);
degR1y2:=finddeg(eq,R1,Y1R1,y2);
degR1y3:=finddeg(eq,R1,Y1R1,y3);
tdegR1:=degR1y1+degR1y2+degR1y3:

degR2y1:=finddeg(eq,R2,Y2R2,y1);
degR2y2:=finddeg(eq,R2,Y2R2,y2);
degR2y3:=finddeg(eq,R2,Y2R2,y3);
tdegR2:=degR2y1+degR2y2+degR2y3:

degR3y1:=finddeg(eq,R3,Y3R3,y1);
degR3y2:=finddeg(eq,R3,Y3R3,y2);
degR3y3:=finddeg(eq,R3,Y3R3,y3);
tdegR3:=degR3y1+degR3y2+degR3y3:

gu:=ptorpolstik(eq,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3,
 ORDER,resh):

#print(`gu is`,gu):

if gu=0 then
 degR1y1:=degR1y1+1:
 degR1y2:=degR1y2+1:
 degR1y3:=degR1y3+1:
 tdegR1:=tdegR1+3:

 degR2y1:=degR2y1+1:
 degR2y2:=degR2y2+1:
 degR2y3:=degR2y3+1:
 tdegR2:=tdegR2+3:

 degR3y1:=degR3y1+1:
 degR3y2:=degR3y2+1:
 degR3y3:=degR3y3+1:
 tdegR3:=tdegR3+3:


gu:=ptorpolstik(eq,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3, ORDER,
resh):
#print(`gu is`,gu):
 if gu=0 then
#  print(`probably ORDER is too low, or S1,S2,S3, are inappropiate, try`):
 # print(`different ones`):
  RETURN(0):
 fi:

gu:=ptorpols(eq,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3,ORDER,N);


if gu=0 then
 degR1y1:=degR1y1+1:
 degR1y2:=degR1y2+1:
 degR1y3:=degR1y3+1:
 tdegR1:=tdegR1+3:

 degR2y1:=degR2y1+1:
 degR2y2:=degR2y2+1:
 degR2y3:=degR2y3+1:
 tdegR2:=tdegR2+3:

 degR3y1:=degR3y1+1:
 degR3y2:=degR3y2+1:
 degR3y3:=degR3y3+1:
 tdegR3:=tdegR3+3:


gu:=ptorpolstik(eq,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3, ORDER,
resh):

 if gu=0 then
#  print(`probably ORDER is too low, or S1,S2,S3, are inappropiate, try`):
# print(`different ones`):
  RETURN(0):
 else
  gu:=ptorpols(eq,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
  degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3,ORDER,N);
   if(gu=0) then
#    print(`probably ORDER is too low, or S1,S2,S3 are inappropiate, try`):
#    print(`different ones`):
    RETURN(0):
   fi:
 fi:
fi:
fi:

if gu=0 then
# print(gu);
# print(`Sorry, try a higher order recurrence`):
RETURN(0):
fi:

if gu=1 then
  gu:=ptorpols(eq,y1,y2,y3,degR1y1,degR1y2,degR1y3,tdegR1,
  degR2y1,degR2y2,degR2y3,tdegR2,degR3y1,degR3y2,degR3y3,tdegR3,ORDER,N)
fi:
ope:=gu[1]:
R_1:=gu[2]*S1:
R_2:=gu[3]*S2:
R_3:=gu[4]*S3:
ope:=normal(ope):
mekh:=denom(ope):
ope:=numer(ope):
ope:=pashet(ope,N):
R_1:=normal(R_1*mekh):
R_2:=normal(R_2*mekh):
R_3:=normal(R_3*mekh):

if ope=0 then
RETURN(0):
else
RETURN(ope,R_1,R_2,R_3):
fi:
end:

#paper:=proc():
#lprint(`Theorem:`):
#lprint(``):
#lprint(`Let G(`,y1,`,`,y2,`,`,y3,`,`,n,`) be`):
#lprint(``):
#print(P*F,``):
#lprint(``):
#lprint(`and a(`,n,`) be its triple integral w.r.t to`,y1,y2,`and`,y3,`.`):
#lprint(``):
#lprint(`Let`,N,`be the forward shift operator in`,n,`.`):
#lprint(``):
#lprint(`The sequence a(`,n,`) satisfies the recurrence`):
#lprint(``):
#print(ope*a(n)=`0.`):
#lprint(``):
#lprint(`Proof: It is routinely verifiable that`):
#lprint(``):
#print(ope*G(y1,y2,y3,n)):
#lprint(``):
#lprint(`= D_`,y1,`(`):
#lprint(``):
#print(R_1*G(y1,y2,y3,n)/P):
#lprint(`          )`):
#lprint(``):
#lprint(`+D_`, y2,`(`):
#lprint(``):
#print(R_2*G(y1,y2,y3,n)/P):
#lprint(``):
#lprint(`+D_`, y3,`(`):
#lprint(``):
#print(R_3*G(y1,y2,y3,n)/P):
#lprint(`          ),`):
#lprint(``):
#lprint(`and the result follows by triple integrating w.r.t to`):
#lprint(y1,`,`,y2,`and`,y3,`.`):
#end:


WZ3:=proc(F,n,x,y,z,resh,N)
local S1,S2,S3,ORDER,k,gu:

for ORDER from 0 to 6 do
for k from 0 to ORDER do
S1:=1/(x*y*z)^k:S2:=1/(x*y*z)^k:S3:=1/(x*y*z)^k:

gu:=WZ3s(1,F,n,x,y,z,ORDER,resh,N,S1,S2,S3):

if gu<>0 then
 RETURN(gu[1],normal(S1*gu[2]),normal(S2*gu[3]),normal(S3(gu[4]))):
fi:

od:
od:

print(`No recurrence of ORDER <=6 found`):
0:
end:


WZ2:=proc(F,n,x,y,resh,N)
local S1,S2,S3,ORDER,k,z,gu:

for ORDER from 0 to 6 do
for k from 0 to ORDER do
S1:=1/(x*y)^k:S2:=1/(x*y)^k:S3:=0:

gu:=WZ3s(1,F,n,x,y,z,ORDER,resh,N,S1,S2,S3):

if gu<>0 then
 RETURN(gu[1],normal(S1*gu[2]),normal(S2*gu[3])):
fi:

od:
od:

print(`No recurrence of ORDER <=6 found`):
0:
end:




degree1:=proc(a,x)
if a=0 then
 RETURN(0):
else
RETURN(degree(a,x)):
fi:
end: