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



print(`Written by Mohamud Mohammed and Doron Zeilberger.`):
print(`First Written: July 22, 2004; tested for Maple 8 and 9. `):
print(`Version of Aug. 10, 2004.  `):
print():
print(`This is ZEILBERGER, A Maple package`):
print(`that implements the new-improved (and simplified) Zeilberger`):
print(`algorithm for single-sum definite hypergeometric summation.`):
print(`It does what zeil does in EKHAD and Zeilberger in the built-in`):
print(`Maple package SumTools[Hypergeometric]. `):
print(`It is one of the two packages`):
print(`that accompany the article: Sharp Upper Bounds for the `):
print(`Orders of the Recurrences Outputted by Zeilberger's algorithm`):
print(`by Mohamud Mohammed and Doron Zeilberger,`):
print(`available from the authors' websites. `):
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 available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name)" .`):
print():
 
 
 
ezra:=proc()
if args=NULL then
print(` ZEILBERGER:  `):
print(` A Maple package for `):
print(` finding recurrences satisfied by definite hypergeomeric sums . `):
print(`Using the Simplified Zeilberger algorithm.`):
print(` For help with a specific procedure, type ezra(procedure_name); `):
print(`The main procedure is  Doron . `):
print():
fi:

if nargs=1 and args[1]=Doron then
print(` Doron(F,k,n,N): Given a hypergeometric term F in the `):
print(` variables k and n, and a symbol N for the shift operator in n:`):
print(` [i.e. Nf(n):=f(n+1)], returns an operator ope(n,N) and`):
print(` a certificate, a rational function R(n,k)`):
print(` such that  G(n,k):=R(n,k)*F(n,k) `):
print(` satisfies ope(N,n)F(k,n)=G(n,k+1)-G(n,k)`):
print(` A hypergeometric term is a product of factorials of the form`):
print(`(an+bk+c)! where a is an integer, b is an integer,`):
print(`and c is symbolic, times a polynomial in (n,k) times x^k for some`):
print(`numeric or symbolic x times y^n for some numeric or symbolic y`):
print(`One can also use binomial(n,k) and RF(a,k) (where RF(a,k))`):
print(`is the rising factorial a(a+1)...(a+k-1) .`):
print(``):
print(` For example, try:`):
print(` Doron(binomial(n,k)^2,k,n,N);`):
print(` Doron(binomial(n,k)^3*2^k*(n^3+k),k,n,N);`):
print(` Doron(1/k!^3/(n-k)!*(n^3+k),k,n,N);`):
print(` Doron(RF(a,k)*RF(n,k),k,n,N);`):
fi:

if nargs=1 and args[1]=preDoron1 then
print(` preDoron1(F,k,n,N): Given a hypergeometric term F in the `):
print(`variables k and n, and a symbol N for the shift operator in n:`):
print(` [i.e. Nf(n):=f(n+1)], returns an operator ope(n,N) and`):
print(` a pre-certificate, a polynomial R(n,k)`):
print(`from which one can construct the certifiate G(n,k)`):
print(` such that ope(N,n)F(k,n)=G(n,k+1)-G(n,k)`):
print(` A hypergeometric term is a product of factorials of the form`):
print(`(an+bk+c)! where a is a NON-negative integer, b is an integer,`):
print(`and c is symbolic, times a polynomial in (n,k) times x^k for some`):
print(`numeric or symbolic x`):
print(`One can also use binomial(n,k) and RF(a,k) (where RF(a,k))`):
print(`is the rising factorial a(a+1)...(a+k-1) .`):
print(``):
print(` For example, try:`):
print(` preDoron1(binomial(n,k)^2,k,n,N);`):
print(` preDoron1(binomial(n,k)^3*2^k*(n^3+k),k,n,N);`):
print(` preDoron1(1/k!^3/(n-k)!*(n^3+k),k,n,N);`):
print(` preDoron1(RF(a,k)*RF(n,k),k,n,N);`):
fi:
end:



ez:=proc():print(`RatioH(Ali,Bli,Cli,Dli,n,L,i) `):
print(`c(Ali,Bli,Cli,Dli,POL,e,n,L)`):
print(` a(Num,Den,z,k,n,L) , b(Num,Den,k,n,L) :`):
print(`Equ1(Num,Den,POL,z,k,n,L,e,X)`):
print(`Z1a(Ali,Bli,Cli,Dli,POL,z,k,n,N)`):
print(`Z1(F,k,n,N),(F=[Num,Den,Pol,z]) `):
print(`Z1L(F,k,n,L,N),(F=[Num,Den,Pol,z]) `):
print(`Z1aL(Ali,Bli,Cli,Dli,POL,z,k,n,L,N)`):
print(`hafokh(F,k,n)`):
print(`preDoron1(F,k,n,N)`):
print(`RecEq(F,k,n,N,low1,high1)`):
end:


RF:=proc(a,n):(a+n-1)!/(a-1)!:end:

rf:=proc(a,n) local i: if n<0 then  0 elif n=0 then  1
else convert([seq(a+i,i=0..n-1)],`*`) fi: end:

RatioH:=proc(Num,Den,n,L,i) local j:
convert([seq( rf(Num[j]+1,i*coeff(Num[j],n,1)),j=1..nops(Num))],`*`)*
convert([seq( rf(
subs(n=n+i,Den[j])+1,(L-i)*coeff(Den[j],n,1)),j=1..nops(Den))],`*`):
end:

c:=proc(Num,Den,POL,e,n,L,w) local i:convert(
[seq(e[i]*subs(n=n+i,POL)*RatioH(Num,Den,n,L,i)*w^i,i=0..L)],`+`):
end:

a:=proc(Num,Den,z,k,n,L) local j,gu:
gu:=z:

for j from 1 to nops(Num) do
 if coeff(Num[j],k,1) > 0 then
  gu:=gu*rf(Num[j]+1,coeff(Num[j],k,1)):
 fi:
od:

for j from 1 to nops(Den) do
 if coeff(Den[j],k,1) < 0 then
  gu:=gu*rf(subs({n=n+L,k=k+1},Den[j])+1,-coeff(Den[j],k,1)):
 fi:
od:
gu:
end:



b:=proc(Num,Den,k,n,L) local j,gu:
gu:=1:

for j from 1 to nops(Num) do
 if coeff(Num[j],k,1) < 0 then
  gu:=gu*rf(subs(k=k+1,Num[j])+1,-coeff(Num[j],k,1)):
 fi:
od:

for j from 1 to nops(Den) do
 if coeff(Den[j],k,1) > 0 then
  gu:=gu*rf(subs({n=n+L},Den[j])+1,coeff(Den[j],k,1)):
 fi:
od:
gu:
end:

Equ1:=proc(Num,Den,POL,z,k,n,L,e,X):
a(Num,Den,z,k,n,L)*X(k+1)-
subs(k=k-1,b(Num,Den,k,n,L))*X(k)-
c(Num,Den,POL,e,n,L):
end:

yafe:=proc(ope,N) local i,ope1:ope1:=expand(ope):
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
end:


FindDeg:=proc(a1,b1,c1,k) local degX,lam,A,B,d:
degX:=degree(c1,k)-max(degree(a1,k),degree(b1,k)):

 if degree(a1,k)=degree(b1,k) and 
coeff(expand(a1),k,degree(a1,k))=coeff(expand(b1),k,degree(b1,k))
  then

  lam:=coeff(expand(b1),k,degree(b1,k)): 

 B:=coeff(expand(b1),k,degree(b1,k)-1):
 A:=coeff(expand(a1),k,degree(a1,k)-1):

    d:=normal((B-A)/lam):

 if type(d,integer) then
   degX:=max(degX+1,d):
 else
  degX:=degX+1:
 fi:
 fi:
degX:
end:

Z1L:=proc(F,k,n,L,N) local e,x,j,degX,X,a1,b1,c1,gu,var1,
var,ope,lu,mu,cert,eq,eqTry,varTry,Num,Den,POL,z,w:
Num:=F[1]: Den:=F[2]: POL:=F[3]: z:=F[4]:w:=F[5]:
ope:=convert([seq(e[j]*N^j,j=0..L)],`+`):
a1:=a(Num,Den,z,k,n,L):
b1:=subs(k=k-1,b(Num,Den,k,n,L)):
c1:=c(Num,Den,POL,e,n,L,w):
degX:=FindDeg(a1,b1,c1,k):
X:=convert([seq(x[j]*k^j,j=0..degX)],`+`):
gu:=expand(a1*subs(k=k+1,X)-b1*X-c1):
var1:={seq(x[j],j=0..degX),seq(e[j],j=0..L)}:
eq:={seq(coeff(gu,k,j),j=0..degree(gu,k))}:

if type(z,symbol) then
eqTry:=subs({z=2, n=5/3},eq):
else
 eqTry:=subs(n=5/3,eq):
fi:
varTry:=solve(eqTry,var1):
if subs(varTry,ope)=0 then
 RETURN(FAIL):
fi:

if type(z,symbol) then
eqTry:=subs({z=3,n=7/5},eq):
 else
eqTry:=subs(n=7/5,eq):
fi:

if subs(varTry,ope)=0 then
 RETURN(FAIL):
fi:

print(`Most likely there is an operator of order`, L, `Please be patinet.`):
var:=solve(eq,var1):


mu:=subs(var,ope),subs(var,X):

if mu[1]=0  then
 RETURN(FAIL):
fi:
if mu[2]=0 then
RETURN(mu):
fi:
lu:=normal(mu[1]/mu[2]):
ope:=yafe(numer(lu),N):
cert:=factor(denom(lu)):
ope,cert:
end:

Z1:=proc(F,k,n,N) local L,gu:
gu:=Z1L(F,k,n,0,N):
for L from 1 while gu=FAIL do
 gu:=Z1L(F,k,n,L,N):
od:
gu:
end:



####converting to the input needed by Z1###
FindFactorial:=proc(Prod1)
local Prod2,i,lu,lu1,lu2:
Prod2:=factor(Prod1):
if type(Prod2,function) and op(0,Prod2)=factorial then
 RETURN(Prod1):
fi:

for i from 1 to nops(Prod2) do
 lu:=op(i,Prod2):
 if type(lu,function) and op(0,lu)=factorial  then
  RETURN(lu):
 fi:

 if type(lu,`^`)  then
   lu1:=op(1,lu):
   lu2:=op(2,lu):

     if type(lu1,function) and op(0,lu1)=factorial then
          if not type(lu2,integer) and lu2>0 then
              RETURN(FAIL):
          else
              RETURN(lu1):
          fi:
     fi: 
 fi:
od:
FAIL:
end:


FindPower:=proc(Prod1,k)
local Prod2,i,lu:
Prod2:=factor(Prod1):
if Prod2=1 then
 RETURN(FAIL):
fi:
if type(Prod2,`^`) and degree(op(2,Prod1),k)=1 then
 RETURN(Prod1):
fi:

for i from 1 to nops(Prod2) do
 lu:=op(i,Prod2):
 if type(lu,`^`) and degree(op(2,lu),k)=1 then
  RETURN(lu):
 fi:
od:
FAIL:
end:



preDoron1:=proc(F,k,n,N) local F1:
F1:=simplify(normal(convert(F,factorial)),power):
F1:=hafokh(F1,k,n):
if F1=FAIL then
 print(`Bad input`):
 RETURN(FAIL):
fi:
Z1(F1,k,n,N):
end:

Doron:=proc(F,k,n,N) local lu, F1,b1,kadima,L,ope,cert,c0:
F1:=normal(convert(F,factorial)):
F1:=hafokh(F1,k,n):
lu:=preDoron1(F,k,n,N):
L:=degree(lu[1],N):
b1:=b(F1[1],F1[2],k,n,L):
kadima:=RatioH(F1[1],F1[2],n,L,0):

ope:=lu[1]:
cert:=normal(lu[2]/kadima*subs(k=k-1,b1)/F1[3]):
c0:=coeff(ope,N,degree(ope,N)):
c0:=coeff(c0,n,degree(c0,n)):
ope:=yafe(ope/c0,N):
cert:=cert/c0:
[ope,cert]:
end:


#RecEq(F,k,n,N,low1,high1): a recurrence equation for
#sum(F(n,k),k=low1..high1) where low1 and high1 have the form
#a+b*n, with b a non0negative integer
RecEq:=proc(F1,k,n,N,low1,high1) local lu,ope,cert,i,i1,gu,F:

F:=convert(F1,factorial):
lu:=Doron(F,k,n,N):

ope:=lu[1]: cert:=lu[2]:

gu:=0:
for i from 0 to degree(ope,N) do
gu:=coeff(ope,N,i)*
(-convert([seq(subs({n=n+i,k=low1+i1},F),i1=0..coeff(low1,n,1)*i-1)],`+`)
+convert([seq(subs({n=n+i,k=high1+i1},F),i1=1..coeff(high1,n,1)*i)],`+`)):
od:

gu:=gu+ subs(k=high1+1,cert)*subs(k=high1+1,F)-
subs(k=low1,cert)*subs(k=low1,F):

gu:=normal(simplify(normal(gu))):
ope,gu:
end:





hafokh:=proc(F,k,n)
local F1,Num,Den,POL,z,F1t,F1b,lu,lu1,lu2,i,w,
PosNum,NegNum,PosDen,NegDen:

F1:=factor(normal(F)):

F1t:=factor(numer(F1)):
F1b:=factor(denom(F1)):

Num:=[]:
lu:=FindFactorial(F1t):
while lu<>FAIL do
lu1:=op(lu):
 if degree(lu1,{n,k})>1 or 
 not type(coeff(lu1,k,1),integer) or
   not type(coeff(lu1,n,1),integer) then
    RETURN(FAIL):
 fi:  
Num:=[op(Num),lu1]:
F1t:=normal(F1t/lu):
lu:=FindFactorial(F1t):
od:

z:=1:
lu:=FindPower(F1t,k):
while lu<>FAIL do
lu1:=op(1,lu):
lu2:=op(2,lu):
if degree(lu2,k)>1 then
  RETURN(FAIL):
else
 z:=z*lu1^coeff(lu2,k,1):
 F1t:=normal(simplify(normal(F1t/lu1^(k*coeff(lu2,k,1))))):
fi:

lu:=FindPower(F1t,k):
od:

w:=1:
lu:=FindPower(F1t,n):
while lu<>FAIL do
lu1:=op(1,lu):
lu2:=op(2,lu):
if degree(lu2,n)>1 then
  RETURN(FAIL):
else
 w:=w*lu1^coeff(lu2,n,1):
 F1t:=normal(simplify(normal(F1t/lu1^(n*coeff(lu2,n,1))))):
fi:

lu:=FindPower(F1t,n):
od:

 

 
POL:=F1t:

Den:=[]:

lu:=FindFactorial(F1b):
while lu<>FAIL do
lu1:=op(lu):
 if degree(lu1,{n,k})>1 or 
 not type(coeff(lu1,k,1),integer) or
   not type(coeff(lu1,n,1),integer) then
    RETURN(FAIL):
 fi:  
Den:=[op(Den),lu1]:


F1b:=normal(F1b/lu):
lu:=FindFactorial(F1b):
od:

lu:=FindPower(F1b,k):
while lu<>FAIL do
lu1:=op(1,lu):
lu2:=op(2,lu):

if degree(lu2,k)>1 then
  RETURN(FAIL):
else
 z:=z/lu1^coeff(lu2,k,1):
 F1b:=normal(simplify(normal(F1b/lu1^(k*coeff(lu2,k,1))))):
fi:

lu:=FindPower(F1b,k):
od:

lu:=FindPower(F1b,n):
while lu<>FAIL do
lu1:=op(1,lu):
lu2:=op(2,lu):

if degree(lu2,k)>1 then
  RETURN(FAIL):
else
 w:=w/lu1^coeff(lu2,n,1):
 F1b:=normal(simplify(normal(F1b/lu1^(n*coeff(lu2,n,1))))):
fi:

lu:=FindPower(F1b,n):
od:

POL:=normal(POL/F1b):


if degree(POL,{n,k})=FAIL then
 RETURN(FAIL):
fi:

PosNum:=[]:
NegNum:=[]:
for i from 1 to nops(Num) do 
if coeff(Num[i],n,1)<0 then
 NegNum:=[op(NegNum),Num[i]]:
else
 PosNum:=[op(PosNum),Num[i]]:
fi:
od:
PosDen:=[]:
NegDen:=[]:
for i from 1 to nops(Den) do 
if coeff(Den[i],n,1)<0 then
 NegDen:=[op(NegDen),Den[i]]:
else
 PosDen:=[op(PosDen),Den[i]]:
fi:
od:
Num:=[op(PosNum),seq(expand((-1)*NegDen[i]-1),i=1..nops(NegDen))]:
Den:=[op(PosDen),seq(expand((-1)*NegNum[i]-1),i=1..nops(NegNum))]:
z:=z*(-1)^(add(coeff(NegDen[i],k,1),i=1..nops(NegDen))
          +add(coeff(NegNum[i],k,1),i=1..nops(NegNum))):

w:=w*(-1)^(add(coeff(NegDen[i],n,1),i=1..nops(NegDen))
          +add(coeff(NegNum[i],n,1),i=1..nops(NegNum))):



[Num,Den,POL,z,w]:

end: