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



print(`Written by Doron Zeilberger.`):
print(`First Written: Oct. 29, 2004; tested for Maple 8 and 9. `):
print(`Version of Oct. 29, 2004.  `):
print():
print(`This is FindHyperGeometric, A Maple package`):
print(`that explains WHY Hypergeometric Closed-Form Evaluations exist`):
print(`(or more generally, why sometimes the recurrence has lower order than it "should").`):
print(`More importantly, it finds (potentially) New (and, in some sense, all) Hypergeometric Series`):
print(`closed-form evaluations of a given hypergeometric summand`):
print(`featuring parameters.`):
print():
print(`It accompanies `):
print(`the article: DECONSTRUCTING the ZEILBERGER algorithm`):
print(`by  Doron Zeilberger,`):
print(`available from the author's website. `):
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 general help, and a list of ALL the functions,`):
print(` type "ezra1();". For specific help type "ezra(procedure_name);" .`):

print():
 

_EnvAllSolutions:=true: 
 

ezra1:=proc()
if args=NULL then
print(`All the procedures are: CleanUp, Doron, DoronMiracles, Equ, Find2F1r, `):
print(`FindHG, FindHGd, FindHGr, FindHGf, FindHGg, FindHGope , FirstMiracle`):
print(`FirstMiracleC, HGtoTerm`):
print():
else
 ezra(args):
fi:
end:

####Sample Famous Summands 
Gauss:=RF(a,k)*RF(b,k)/RF(c,k)/k!:
Kummer:=RF(a,k)*RF(b,k)/RF(1+a-b,k)/k!*(-1)^k:
GaussHalf:=RF(2*a,k)*RF(2*b,k)/RF(a+b+1/2,k)/k!*(1/2)^k:
SaalschutzG:=r->RF(a,k)*RF(b,k)*RF(c,k)/k!/RF(d,k)/RF(a+b+c-d+r,k)*x^k:
Saalschutz:=RF(a,k)*RF(b,k)*RF(c,k)/k!/RF(d,k)/RF(a+b+c-d+1,k):
Dixon:=RF(a,k)*RF(b,k)*RF(c,k)/k!/RF(1+a-b,k)/RF(1+a-c,k):
DixonG:=RF(-n,k)*RF(a,k)*RF(c,k)/k!/RF(b+n,k)/RF(d,k)*x^k:
Dougall:=subs(f=1+2*a-b-c-d-e,(2*k+a)*RF(a,k)*RF(b,k)*RF(c,k)*RF(d,k)*RF(e,k)*RF(f,k)/
RF(1+a-b,k)/RF(1+a-c,k)/RF(1+a-d,k)/RF(1+a-e,k)/RF(1+a-f,k)/k!):
Apery:=binomial(n,k)^2*binomial(n+k,k)^2:
####End Famous Summands

ezra:=proc()
if args=NULL then
print(` FindHyperGeometric:  `):
print(` A Maple package for explaining (from the viewpoint of Zeilberger's algorithm)`):
print(` Exact Evaluations  of hypergeometric series, and for discovering`):
print(` (and of course, proving) them from scratch.  `):
print(``):
print(` For help with a specific procedure, type ezra(procedure_name); `):
print(`The main procedures are: DoronMiracles`):
print(`FindHG  , FindHGd, FindHGr , FirstMiracleC `):
print():
fi:

if nargs=1 and args[1]=CleanUp then
print(`CleanUp(kv,n): cleans up equivalent HGs`):
fi:

if nargs=1 and args[1]=Equ then
print(`Equ(HG1,HG2),n: are the Hypergeometric terms`):
print(`HG1 and HG2 (phrased in terms of n) equivalent?`):
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]=DoronMiracles then
print(` DoronMiracles(F,k,n,L,N): a verbose form of Doron (q.v) where the`):
print(`order is known, explaining all the miracles, if any.`):
print(`Gauss, GaussHalf, Kummer, Dixon, Saalschutz, Dougall and Apery, are built-in.`):
print(`To understand the miracles behind Guass's 2F1(a,b;c;1), Guass's x=1/2, Kummer's`):
print(`2F1(a,n;1+a-n;-1), Pfaff-Saalschutz's 3F2(a,b,c;d,a+b+c-d+1;1)`);
print(` Dixon's 3F2(a,b,c;1+a-b,1+a-c;1) and Dougall's 7F6, type, respectively`):
print(` DoronMiracles(Gauss,k,a,1,A);`):
print(` DoronMiracles(GaussHalf,k,a,1,A);`):
print(` DoronMiracles(Kummer,k,b,1,B);`):
print(` DoronMiracles(Saalschutz,k,a,1,A);`):
print(` DoronMiracles(Dixon,k,b,1,B);`):
print(` DoronMiracles(Dougall,k,b,1,B);`):
print(`To understand the miracle behind Apery's recurrence, type`):
print(` DoronMiracles(Apery,k,n,2,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 certificate 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:




if nargs=1 and args[1]=Find2F1r then
print(`Find2F1r(K,n,R): compiles a table such that`):
print(`T[a,b,c] is the set of successful closed-form`):
print(`evaluations of 2F1([-a*n,b*n+beta],[c*n+gamma],z)`):
print(`with 1<=a<=K, -K<= b <= K, -K <= c <= K c<>{a,b}`):
print(`and not consequences of previous ones`):
print(`R is the symbol for rising-factorial`):
print(`try Find2F1r(1,n,R)`):
fi:

if nargs=1 and args[1]=FindHG then
print(`FindHG(HG,n,L,Pars): Given a hypergeometric term`):
print(`HG=[Num,Den,z] template, phrased in terms`):
print(`of n and the parameters Pars, finds the specializations that`):
print(`yield recurrences of order L, thanks to a miracle of the`):
print(`third kind. For example, try:`):
print(`FindHG([[-n,-n+a],[2*n+b],z],n,1,{a,b,z});`):
fi:

if nargs=1 and args[1]=FindHGr then
print(`FindHGr(HG,n,Pars,R): Given a hypergeometric term`):
print(`HG=[Num,Den,z] template, phrased in terms`):
print(`of n and the parameters Pars, finds the specializations that`):
print(`yields Closed-Form evaluations, in terms of rising factorials`):
print(`thanks to a miracle of the`):
print(`third kind. R is the symbol for rising factorial`):
print(`i.e. R(a,n):=a(a+1) ... (a+n-1) `):
print(`For example, try:`):
print(`FindHGr([[-n,-n+a],[2*n+b],z],n,{a,b,z},R);`):
fi:

if nargs=1 and args[1]=FindHGf then
print(`FindHGf(HG,n,Pars): Given a hypergeometric term`):
print(`HG=[Num,Den,z] template, phrased in terms`):
print(`of n and the parameters Pars, finds the specializations that`):
print(`yields Closed-Form evaluations (n terms of factorials),`):
print(` thanks to a miracle of the`):
print(`third kind.`):
print(`For example, try:`):
print(`FindHGf([[-n,-n+a],[2*n+b],z],n,{a,b,z});`):
fi:

if nargs=1 and args[1]=FindHGg then
print(`FindHGg(HG,n,Pars): Given a hypergeometric term`):
print(`HG=[Num,Den,z] template, phrased in terms`):
print(`of n and the parameters Pars, finds the specializations that`):
print(`yields Closed-Form evaluations (n terms of GAMMA),`):
print(` thanks to a miracle of the`):
print(`third kind. R is the symbol for rising factorial`):
print(`i.e. R(a,n):=a(a+1) ... (a+n-1) `):
print(`For example, try:`):
print(`FindHGg([[-n,-n+a],[2*n+b],z],n,{a,b,z});`):
fi:

if nargs=1 and args[1]=FindHGope then
print(`FindHGope(HG,n,L,Pars,N): Given a hypergeometric term`):
print(`HG=[Num,Den,z] template, phrased in terms`):
print(`of n and the parameters Pars, finds the specializations that`):
print(`yield recurrences of order L, thanks to a miracle of the`):
print(`[followed by the operator where N is the shift operator in n]`):
print(`third kind. For example, try:`):
print(`FindHGope([[-n,-n+a],[2*n+b],z],n,1,{a,b,z},N);`):
fi:

if nargs=1 and args[1]=FindHGd then
print(`FindHGd(F,k,n,L,Pars): given a  hypergeometric term (Directly in terms of`):
print(` factorials or rising factorials RF(a,k)) `):
print(`with parameters Pars, finds for what specializations ot these parameters,`):
print(`sum(F,k) has a Z-certificate with  a recurrence of order L `):
print(`For example, for discovering the Pfaff-Saalschutz identity from scratch, try:`):
print(`FindHGd(RF(a,k)*RF(b,k)*RF(c,k)/k!/RF(d,k)/RF(a+e,k),k,a,1,{b,c,d,e});`):
fi:

if nargs=1 and args[1]=FirstMiracle then
print(`FirstMiracle(F1,k,n,L): Does the first miracle happen?`):
print(`for the summand F1 and order L?`):
print(`e.g. try:`):
print(`FirstMiracle(binomial(n,k)^2*binomial(n+k,k)^2,k,n,2);`):
fi:

if nargs=1 and args[1]=FirstMiracleC then
print(`FirstMiracleC(F1,k,n,L,Pars): Given a `):
print(`hypergeometric term F1`):
print(`featuring parameters in Pars, find the specializations`):
print(`that make the first miracle happen`):
print(`for order L?`):
print(`e.g. try:`):
print(`FirstMiracleC(RF(-n,k)*RF(a,k)*`):
print(`RF(b,k)/k!/RF(c,k)/RF(-n+d,k)*x^k,k,n,1,{d,x});`):
print(``):
print(`FirstMiracleC(SaalchutzG(2)*x^k,k,n,1,{x});`);
print(`FirstMiracleC(SaalchutzG(3)*x^k,k,n,1,{x});`);



fi:

if nargs=1 and args[1]=HGtoTerm then
print(`HGtoTerm(HG,k): converts F(HG,k) (where HG=[Num,Den,z]) to`):
print(` closed-form`):
print(`For example try HGtoTerm([[-n,a],[b],z],k)`):
fi:

end:


####Start from ZEILBERGER
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))}:


 eqTry:=subs(n=5/3,eq):

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

eqTry:=subs(n=7/5,eq):

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

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:



#####End from ZEILBERGER

#Begin stuff for FindHG

#FindHGd(F,k,n,L,Pars): given a  hypergeometric term with
#parameters Pars, finds for what specializations
#sum(F,k,n) satisfies a recurrence of order L
#For example, try FindHGd(RF(-n,k)*RF(a,k)/k!/RF(b,k)*z^k,k,n,1,{a,b,z});

FindHGd:=proc(F1,k,n,L,Pars) local e,x,j,degX,X,a1,b1,c1,gu,var1,
eq,Num,Den,POL,z,w,Mat1,i1,j1,F,Eq,i2:

F:=hafokh(F1,k,n):

Num:=F[1]: Den:=F[2]: POL:=F[3]: z:=F[4]:w:=F[5]:
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 nops(eq)>=nops(var1) then

 Eq:={}:
    for i2 from nops(var1) to nops(eq) do
  Mat1:=
 [seq([seq(coeff(eq[i1],var1[j1],1),j1=1..nops(var1))],i1=1..nops(var1)-1)
 ,[seq(coeff(eq[i2],var1[j1],1),j1=1..nops(var1))]]:

  gu:=expand(linalg[det](Mat1)):
  Eq:=Eq union {seq(coeff(gu,n,i1),i1=0..degree(gu,n))}:
 
 od:
fi:
{solve(Eq,Pars)}:
end:


####Verbose

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

 if degree(a1,k)=degree(b1,k) then

print(`The leading degrees match`):


ka:=expand(coeff(expand(a1),k,degree(a1,k))-coeff(expand(b1),k,degree(b1,k))):
print(`the difference between the leading coefficients is:`,ka):

if ka<>0 then
print(`A miracle of the First kind didn't happen.`):
 else

   print(`A miracle of the first kind happened, hence the degree of X(k) is one higher,namley`, degX+1):
  
    if degX+1<0 then
      print(`Since it is negative, there is no hope, unless a miracle of the second kind would happen.`):
   fi:
  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):


print(` Possible Extra degree(beyond degX+1) is ` ,d-degX-1):

if type(d,integer) then

 if d>degX+1 then
   print(`A miracle of the second kind happened!`):
   degX:=max(degX+1,d):
   print(`The degree of X(k) is yet higher, namely`, d):
 else
    print(`A miracle of the second kind did not happen, since the alternative degree is not higher.`):
  degX:=degX+1:
 fi:



else
 degX:=degX+1:
   print(`A miracle of the second kind did not happen, since the alternative degree is symbolic.`):

fi:
fi:
fi:

degX:
end:

DoronMiracles:=proc(F1,k,n,L,N) local e,x,j,degX,X,a1,b1,c1,gu,var1,
var,ope,lu,mu,cert,eq,Num,Den,POL,z,w,F,Lg:
print(`The summand is`, F1):
F:=hafokh(convert(F1,factorial),k,n):
Num:=F[1]: Den:=F[2]: POL:=F[3]: z:=F[4]:w:=F[5]:

Lg:=GenericOrder(Num,Den,k):

print(`The guaranteed (generic) upper bound for the Order is: `,Lg):

if L<Lg then
print(`Since the order that you are after`, L, `is less than that,`):
print(`in order for a recurrence of that order to happen, some`):
print(`miracle(s) are bound to happen.`):
else
 print(`This can be achieved without miracles.`):
fi:

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:=FindDegV(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))}:



print(`The number of equations is`,nops(eq)):
print(`The number of unknowns is`,nops(var1)):
var:=solve(eq,var1):


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

if mu[1]=0  then
 print(`Sorry, there are no miracles, and there is no (Z-)recurrence of order`, L):
 RETURN(FAIL):
fi:
if mu[2]=0 then
RETURN(mu):
fi:

if nops(eq)>=nops(var1) then
 print(`A miracle of the third kind happened`):
fi:

lu:=normal(mu[1]/mu[2]):
ope:=yafe(numer(lu),N):
cert:=factor(denom(lu)):
ope,cert:
end:


###


#HGtoTerm(HG,k): converts F(HG,k) (where HG=[Num,Den,z]) to closed-form
#For example try GHtoTerm([[-n,a],[b],z],k)
HGtoTerm:=proc(HG,k)
local i:
convert([seq(RF(HG[1][i],k),i=1..nops(HG[1]))],`*`)/
convert([seq(RF(HG[2][i],k),i=1..nops(HG[2]))],`*`)/k!*HG[3]^k:
end:


#FindHG(HG,n,L,Pars): Given a hypergeometric term
#HG=[Num,Den,z] template, phrased in terms
#of n and the parameters Pars, finds the specializations that
#yield recurrences of order L, thanks to a miracle of the
#third kind. For example, try:
#FindHG([[-n,-n+a],[2*n+b],z],n,1,{a,b,z});
FindHG:=proc(HG,n,L,Pars) local gu,k,i:
gu:=FindHGd(HGtoTerm(HG,k),k,n,L,Pars):
gu:={seq(subs(gu[i],HG),i=1..nops(gu))}:
CleanUp(gu,n):

end:


Equ1:=proc(HG1,HG2,n,i)
local m,m1,HG1a:
m1:=solve(subs(n=m,HG2[1][1])=HG1[1][i],n):

if {m1}={} then
  RETURN(false):
fi:

HG1a:=expand(subs(m=n,subs(n=m1,HG1))):

if convert(HG1a[1],set)=convert(HG2[1],set) and
  convert(HG1a[2],set)=convert(HG2[2],set)  then
   true:
else
   false:
fi:


end:

Equ:=proc(HG1,HG2,n) local i:

if HG1[3]<>HG2[3] then
 RETURN(false):
fi:

if nops(HG1[1])<>nops(HG2[1]) or
   nops(HG1[2])<>nops(HG2[2]) then
   RETURN(false):
fi:

for i from 1 to nops(HG1[1]) do
 if Equ1(HG1,HG2,n,i) then
   RETURN(true):
 fi:
od:
false:
end:




EquS:=proc(HG1,kv,n)
local i:
 for i from 1 to nops(kv) do
   if Equ(HG1,kv[i],n) then
      RETURN(true):
   fi:
 od:
false:
end:

#CleanUp(kv,n): cleans up equivalent HGs
CleanUp:=proc(kv,n)
local gu,i,kv1:



kv1:={}:

for i from 1 to nops(kv) do
 if not IsBad(kv[i]) then
   kv1:=kv1 union {kv[i]}:
 fi:
od:

if kv1={}  then
 RETURN({}):
fi:


gu:={kv1[1]}:

for i from 2 to nops(kv1) do
 if not EquS(kv1[i],gu,n) then
     gu:=gu union {kv1[i]}:
 fi:
od:
gu:


end:

#FindHGope(HG,n,L,Pars,N): Given a hypergeometric term
#HG=[Num,Den,z] template, phrased in terms
#of n and the parameters Pars, finds the specializations that
#yield recurrences of order L, thanks to a miracle of the
#third kind. For example, try:
#FindHGope([[-n,-n+a],[2*n+b],z],n,1,{a,b,z},N);
FindHGope:=proc(HG,n,L,Pars,N) local gu,k,i:
gu:=FindHGd(HGtoTerm(HG,k),k,n,L,Pars):
gu:={seq(subs(gu[i],HG),i=1..nops(gu))}:
gu:=CleanUp(gu,n):
{seq([gu[i],Doron(HGtoTerm(gu[i],k),k,n,N)[1]],i=1..nops(gu))}:
end:


#PtorOpe(ope,n,N,R): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of rising factorial R(n,k) (R is a symbol)
#For example try PtorOpe((n+1)*N-1,n,N,R);
PtorOpe:=proc(ope,n,N,R)
local rat,c,t1,b1,t1c,b1c,i:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

c:=t1c/b1c:

t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

c^n* convert([seq(R(-t1[i],n),i=1..nops(t1))],`*`)/
convert([seq(R(-b1[i],n),i=1..nops(b1))],`*`):

end:

#PtorOpeF(ope,n,N): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of factorials (R is a symbol)
#For example try PtorOpeF((n+1)*N-1,n,N,R);
PtorOpeF:=proc(ope,n,N)
local rat,c,t1,b1,t1c,b1c,i:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

c:=t1c/b1c:

t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

for i from 1 to nops(t1) do
 if type(t1[i],integer) and t1[i]>=0 then
   RETURN(FAIL):
 fi:
od:

for i from 1 to nops(b1) do
 if type(b1[i],integer) and b1[i]>=0 then
   RETURN(FAIL):
 fi:
od:

c^n* 
normal(convert([seq((-t1[i]+n-1)!/(-t1[i]-1)!,i=1..nops(t1))],`*`)/
convert([seq((-b1[i]+n-1)!/(-b1[i]-1)!,i=1..nops(b1))],`*`)):

end:



FindHGr:=proc(HG,n,Pars,R) local gu,N,i,mu,lu1,lu1F,ope,cert,coe,k,mu1:
gu:=FindHGope(HG,n,1,Pars,N):
mu:={}:
for i from 1 to nops(gu) do
lu1:=PtorOpe(gu[i][2],n,N,R):
lu1F:=PtorOpeF(gu[i][2],n,N):

if lu1F<>FAIL then
ope:=Doron(normal(HGtoTerm(gu[i][1],k)/lu1F),k,n,N):
cert:=ope[2]:
ope:=expand(ope[1]):
coe:=coeff(ope,N,1):

if coe<>0 then
ope:=factor(normal(ope/coe)):
cert:=normal(cert/coe):

ope:=factor(coeff(ope,N,1))*N+factor(coeff(ope,N,0)):

if ope<>N-1 then
 ERROR(`Something is wrong`):
fi:

mu1:=[gu[i][1],lu1,cert]:


 mu:=mu union {mu1}:
fi:
fi:
od:
mu:
end:


FindHGf:=proc(HG,n,Pars) local gu,N,i,mu,lu1,ope,cert,coe,k,mu1:
gu:=FindHGope(HG,n,1,Pars,N):
mu:={}:
for i from 1 to nops(gu) do
lu1:=PtorOpeF(gu[i][2],n,N):

if lu1<>FAIL then
ope:=Doron(normal(HGtoTerm(gu[i][1],k)/lu1),k,n,N):
cert:=ope[2]:
ope:=expand(ope[1]):
coe:=coeff(ope,N,1):
ope:=factor(normal(ope/coe)):
cert:=normal(cert/coe):

ope:=factor(coeff(ope,N,1))*N+factor(coeff(ope,N,0)):

if ope<>N-1 then
 ERROR(`Something is wrong`):
fi:

mu1:=[gu[i][1],lu1,cert]:


 mu:=mu union {mu1}:
fi:
od:
mu:
end:



IsBad:=proc(HG)
local i:

for i from 1 to nops(HG[1]) do
 if type(HG[1][i],integer) and HG[1][i]<0 then
   RETURN(true):
 fi:
od:

for i from 1 to nops(HG[2]) do
 if type(HG[2][i],integer) and HG[2][i]<0 then
   RETURN(true):
 fi:
od:

false:
end:


FindHGg:=proc(HG,n,Pars) local gu,N,i,mu,lu1,ope,cert,coe,k,mu1:
gu:=FindHGope(HG,n,1,Pars,N):
mu:={}:
for i from 1 to nops(gu) do
lu1:=PtorOpeF(gu[i][2],n,N):

if lu1<>FAIL then
ope:=Doron(normal(HGtoTerm(gu[i][1],k)/lu1),k,n,N):
cert:=ope[2]:
ope:=expand(ope[1]):
coe:=coeff(ope,N,1):
ope:=factor(normal(ope/coe)):
cert:=normal(cert/coe):

ope:=factor(coeff(ope,N,1))*N+factor(coeff(ope,N,0)):

if ope<>N-1 then
 ERROR(`Something is wrong`):
fi:

mu1:=[gu[i][1],convert(lu1,GAMMA),cert]:


 mu:=mu union {mu1}:
fi:
od:
mu:
end:




#Find2F1r(K,n,R): compiles a table such that
#T[a,b,c] is the set of successful closed-form
#evaluations of 2F1([-a*n,b*n+beta],[c*n+gamma],z)
#with 1<=a<=K, -K<= b <= K, -K <= c <= K c<>{a,b}
#and not consequences of previous ones
#R is the symbol for rising-factorial
#try Find2F1r(1,n,R)
Find2F1r:=proc(K,n,R)
local T,gu,mu,b1,c1,z,mu1,a,b,c,gu1,i,ku:
ku:={}:
gu:={}:
gu1:={}:
for a from 1 to K do
for b from -K to K do
 for c from -K to K do

  if c<>a and c<>b then
   mu1:=FindHGr([[-a*n,b*n+b1],[c*n+c1],z] ,n,{b1,c1,z},R):
    mu:={}:
     for  i from 1 to nops(mu1) do
      if not EquS(mu1[i][1],gu1,n) then
        mu:=mu union {mu1[i]}:
        gu1:=gu1 union {mu1[i][1]}:
        gu:=gu union {mu1[i]}:
      fi:
     od:
  if mu<>{} then
      ku:=ku union {[-a,b,c]}:
       T[-a,b,c]:=mu:
  fi:
  fi:
  
od:od:od:  
ku,op(T):
end:




#GenericOrder(Num,Den,k): the generic order 
GenericOrder:=proc(Num,Den,k)
local Sa,Sb,Sc,Sd,i:

Sa:=0: Sb:=0: Sc:=0: Sd:=0:

for i from 1 to nops(Num) do
 if coeff(Num[i],k,1)>=0 then
   Sa:=Sa+coeff(Num[i],k,1):
 else
   Sb:=Sb-coeff(Num[i],k,1):
 fi:
od:

for i from 1 to nops(Den) do
 if coeff(Den[i],k,1)>=0 then
   Sc:=Sc+coeff(Den[i],k,1):
 else
   Sd:=Sd-coeff(Den[i],k,1):
 fi:
od:

max(Sa+Sd,Sb+Sc):
end:


FirstMiracle1:=proc(a1,b1,k) local ka:

 if degree(a1,k)<>degree(b1,k) then
  RETURN(false):
 fi:

ka:=expand(coeff(expand(a1),k,degree(a1,k))-coeff(expand(b1),k,degree(b1,k))):

if ka=0 then
 RETURN(true):
else
 RETURN(false):
fi:

end:

#FirstMiracle(F1,k,n,L): Does the first miracle happen?
#for the summand F1 and order L?
#e.g. try:
#FirstMiracle(binomial(n,k)^2*binomial(n+k,k)^2,k,n,2);
FirstMiracle:=proc(F1,k,n,L) local a1,b1,Num,Den,z,F,Lg:
F:=hafokh(convert(F1,factorial),k,n):
Num:=F[1]: Den:=F[2]: z:=F[4]:

Lg:=GenericOrder(Num,Den,k):

if L>=Lg then
 print(`You don't need miracles`):
  RETURN(FAIL):
fi:

a1:=a(Num,Den,z,k,n,L):
b1:=subs(k=k-1,b(Num,Den,k,n,L)):
FirstMiracle1(a1,b1,k):

end:



FirstMiracle1C:=proc(a1,b1,k,Pars) local ka:

 if degree(a1,k)<>degree(b1,k) then
  RETURN(FAIL):
 fi:

ka:=expand(coeff(expand(a1),k,degree(a1,k))-coeff(expand(b1),k,degree(b1,k))):

ka:={solve(ka,Pars)}:

end:


#FirstMiracleC(F1,k,n,L,Pars): Given a hypergeometric term F1
#featuring parameters in Pars, find the specializations
#that make the first miracle happen
#for order L?
#e.g. try:
#FirstMiracle(RF(-n,k)*RF(a,k)*RF(b,k)/k!/RF(c,k)/RF(-n+d,k)*x^k,k,n,1.{d,x});
FirstMiracleC:=proc(F1,k,n,L,Pars) local a1,b1,Num,Den,z,F,Lg,ka,i:
F:=hafokh(convert(F1,factorial),k,n):
Num:=F[1]: Den:=F[2]: z:=F[4]:

Lg:=GenericOrder(Num,Den,k):

if L>=Lg then
 print(`You don't need miracles`):
  RETURN(FAIL):
fi:

a1:=a(Num,Den,z,k,n,L):
b1:=subs(k=k-1,b(Num,Den,k,n,L)):
ka:=FirstMiracle1C(a1,b1,k,Pars):

if ka=FAIL then
 RETURN({}):
fi:

{seq([ka[i],subs(ka[i],F1)],i=1..nops(ka))}:
end:







PfaffT1:=proc(HG):[HG[1],[HG[1][1]+HG[1][2]+1-HG[2][1]],1-HG[3]]:end: