#######################################################################
## twoFone Save this file as twoFone. To use it                       #
## stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read twoFone:             <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 Version: Oct. 29, 2004; tested for Maple 8 and 9. `):
print(`This Version: Oct. 29, 2004.  `):
print():
print(`This is twoFone, A Maple package`):
print(`that finds (in some sense, all) 2F1 Hypergeometric Series`):
print(`closed-form evaluations featuring parameters.`):
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 a list of all functions, type ezra1(); `):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name)" .`):
print():
 

_EnvAllSolutions:=true: 
 

ezra1:=proc()
if args=NULL then
print(`All the procedures are: AZd,Book2F1,BOOK, CandSet,CleanUp, Doron, DoronD, DoronV, Equ, Euler, `):
print(`Find2F1, Find2F1R, Find2F1Ra, Find2F1wz `):
print(`FindHG, FindHGc, FindHGr, FindHGr1,FindHGf, FindHGg, FindHGope, FindHG1 , FirstMiracle`):
print(`FirstMiracleC, HGtoTerm, HypToInt1, HypToInt`):
print(`Implications, ImplicationsL, IsMultiple , Khaverim1, `):
print(`PfaffT1, PfaffT2, Quad1, Quad2, Recip, Yafe`):
print():
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then
print(` twoFone :   `):
print(` A Maple package for `):
print(` finding Exact Evaluations of  `):
print(` Hypergeometric Series Closed-Form Evaluations`):
print(` For help with a specific procedure, type ezra(procedure_name); `):
print(`The main procedures are:  BOOK, Doron, Find2F1wz, FindHG`):
print(`  FindHGr, FindHGc. `):
print():
fi:

if nargs=1 and args[1]=AZd then
print(`AZd(A,y,n,N): performs the Almkvist-Zeilberger algorithm; it finds a`):
print(`recurrence,in the discrete variable n, phrased in terms of the shift operator N`):
print(`for int(A(y,n),y), where  A(y,n) is hypergeometric-hyperexponential in (y,n) `):
print(`For example, try: AZd(exp(-x)*x^n/n!,x,n,N);`):
fi:

if nargs=1 and args[1]=Book2F1 then
print(`Book2F1(K,n,R,L): writes a book about exact evaluations of 2F1's of the form`):
print(`2F1(-a*n,b*n+b1;c*n+c1;z) with  1<=a<=K, -K<=b,c<=K, complete with proofs (certificates)`):
print(`It finds ALL non-trivial WZ-able identities in the range that are NOT consequences`):
print(`of the classical identities of Gauss, and Kummer, and omits consequences of previous`):
print(`identities by Euler and Pfaff's transformations (up to L iterations)`):
print(`R is the symbol for the rising factorial, i.e.`):
print(` R(a,k):=a(a+1)(a+2)...(a+k-1)`):
print(``):
print(`For example, try: Book2F1(2,n,R,3);`):
fi:

if nargs=1 and args[1]=BOOK then
print(`BOOK(K,K1,n,R,L): writes a book about exact evaluations of 2F1's of the form`):
print(`2F1(-a*n,b*n+b1;c*n+c1;z) with  1<=a<=K1, -K<=b,c<=K, complete with proofs (certificates)`):
print(`It finds ALL non-trivial WZ-able identities in the range that are NOT consequences`):
print(`of the classical identities of Gauss, and Kummer, and omits consequences of previous`):
print(`identities by Euler and Pfaff's transformations (up to L iterations)`):
print(`R is the symbol for the rising factorial, i.e.`):
print(` R(a,k):=a(a+1)(a+2)...(a+k-1)`):
print(``):
print(`For example, try: BOOK(2,1,n,R,3);`):
fi:

if nargs=1 and args[1]=CandSet then
print(`CandSet(K,L): a;; the triples of integers between`):
print(`-K and K that are inequivalent to each other up to`):
print(`L degrees of separation`):
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): performs the simplified Zeilberger algorithm`):
print(`(for PROPER hypergeometric summation) as described in Mohammed and Zeilberger's`):
print(`article: "Sharp Upper Bounds for the Recurrences outputted by the `):
print(`Zeilberger and q-Zeilberger algorithms."`):
print(``):
print(`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]=DoronD then
print(`DoronD(F,k,n,N): Does Doron with F given in hypergeometric notation [Numer,Den,z]`):
print(`For example, try: DoronD([[-n,a],b,1],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 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]=Euler then
print(`Euler(HG): Applies (2.2.7) in AAR`):
fi:


if nargs=1 and args[1]=Find2F1 then
print(`Find2F1(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 trivial consequences of previous ones`):
print(`R is the symbol for rising-factorial`):
print(`try Find2F1(1,n,R)`):
fi:

if nargs=1 and args[1]=Find2F1R then
print(`Find2F1R(K,n,R,L): 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 (by Euler, PfaffT1, PfaffT2, Recip, Symm`):
print(` up to L iterations`): 
print(`R is the symbol for rising-factorial`):
print(`try Find2F1R(2,n,R,5)`):
fi:

if nargs=1 and args[1]=Find2F1Ra then
print(`Find2F1Ra(K,K1,n,R,L): 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<=K1, -K<= b <= K, -K <= c <= K c<>{a,b}`):
print(`and not consequences of previous ones (by Euler, PfaffT1, PfaffT2, Recip, Symm`):
print(` up to L iterations`): 
print(`R is the symbol for rising-factorial`):
print(`try Find2F1Ra(2,2,n,R,5)`):
fi:

if nargs=1 and args[1]=Find2F1wz then
print(`Find2F1wz(K,K1,n,L,y): Finds all discrete-cont.`):
print(`WZ pairs [F(n,y),G(n,y)] (i.e. that satisfy`):
print(`F(n+1,y)-F(n,y)=diff(G(n,y),y) `):
print(`that come from the hypergeometric closed-form evaluations that come`):
print(` from`):
print(`Find2F1(K,K1,n,R,L) (q.v.`):
print(`try Find2F1wz(2,2,n,5,y)`):
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]=FindHGc then
print(`FindHGc(HG,n,Pars,R,y): Discrete-Continuous WZ-pairs implied by FindHGr`):
print(`Try: FindHGc([[-n,-n+a],[2*n+b],z],n,{a,b,z},y);`):
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(`The output is a set of lists, where every list is of the form`):
print(`[Succesful HG, RHS (in terms of R's), certificate]`):
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]=FindHGr1 then
print(`FindHGr1(tri,n,R): Given a triplet, tri1:[a,b,c], say`):
print(`Here n is the symbol for the n variable, and R is short for rising factorial`):
print(`Finds all the closed form evaluations (thanks to WZ) of the form`):
print(`2F1(a*n1,b*n+b1;c*n+c1;z) for some b1,c1,z`):
print(`For example, try: FindHGr1([-1,2,0],n,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. R is the symbol for rising factorial`):
print(`i.e. R(a,n):=a(a+1) ... (a+n-1) `):
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]=FindHG1 then
print(`FindHG1(F,k,n,L,Pars): given a  hypergeometric term with`):
print(`parameters Pars, finds for what specializations`):
print(`sum(F,k) satisfies a recurrence of order L`):
print(`For example, try:`):
print(` FindHG1(RF(-n,k)*RF(a,k)/k!/RF(b,k)*z^k,k,n,1,{a,b,z});`):
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});`):
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:

if nargs=1 and args[1]=Implications then
print(`Implications(HG): the set of implications of HG by Recip, Euler, PfaffT1,and PfaffT2`):
fi:

if nargs=1 and args[1]=ImplicationsL then
print(`ImplicationsL(HG): the set of implications of HG by Recip, Euler, PfaffT1,and PfaffT2`):
print(` after L iterations`):
fi:

if nargs=1 and args[1]=HypToInt then
print(`HypToInt(Hyp,Rhs,y): converting a 2F1 Hyp to integral form`):
print(`with Rhs given in terms of R(a,k) denoting rising factorial`):
print(`For example, to convert Gauss's 2f1(a.b,c,1) to an integral, try:`):
print(`HypToInt([[a,b],[c],1],(c-1)!*(c-a-b-1)!/(c-a-1)!/(c-b-1)!,y);`):
fi:

if nargs=1 and args[1]=HypToInt1 then
print(`HypToInt1(Hyp,y): converting a 2F1 to integral form`):
fi:


if nargs=1 and args[1]=IsMultiple then
print(`IsMultiple(HG,n): Is it the case n=r*n for some r`):
fi:

if nargs=1 and args[1]=Khaverim1 then
print(`Khaverim1(tri1): given a triple tri1=[A,B,C], finds`):
print(`all the immediately`):
print(`implied triples by Neg. Symmm. Recipr. Pfaff2, Euler and Pfaff1`):
fi:

if nargs=1 and args[1]=PfaffT1 then
print(`PfaffT1(HG):Applies (2.3.14) in AAR`):
fi:

if nargs=1 and args[1]=PfaffT2 then
print(`PfaffT1(HG):Applies (2.2.6) in AAR`):
fi:

if nargs=1 and args[1]=Quad1 then
print(`Quad1(HG): Applies (3.1.2) in AAR`):
fi:

if nargs=1 and args[1]=Quad2 then
print(`Quad2(HG): Applies (3.1.7) in AAR`):
fi:

if nargs=1 and args[1]=Recip then
print(`Recip(HG): the reciprocal 2F1 of HG`):
fi:

if nargs=1 and args[1]=Yafe then
print(`Yafe(i,HG): writes the theorem given by HG nicely`):
fi:

end:

Saa:=r->RF(-n,k)*RF(a,k)*RF(b,k)/k!/RF(c,k)/RF(-n+r+a+b-c,k):
Saag:=RF(-n,k)*RF(a,k)*RF(b,k)/k!/RF(c,k)/RF(-n+d,k)*x^k:
Dix:=RF(-n,k)*RF(a,k)*RF(c,k)/k!/RF(1+a+n,k)/RF(1+a-c,k);
Dixg:=RF(-n,k)*RF(a,k)*RF(c,k)/k!/RF(b+n,k)/RF(d,k)*x^k;


####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

#FindHG1(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 FindHG1(RF(-n,k)*RF(a,k)/k!/RF(b,k)*z^k,k,n,1,{a,b,z});

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





#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:=FindHG1(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:=FindHG1(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:




#Find2F1(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 Find2F(1,n,R)
Find2F1:=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(HG):Applies (2.3.14) in AAR
PfaffT1:=proc(HG):
if not (nops(HG)=3 and nops(HG[1])=2 and nops(HG[2])=1) then
 print(`Not a 2F1`):
 RETURN(FAIL):
fi:
[HG[1],[HG[1][1]+HG[1][2]+1-HG[2][1]],1-HG[3]]:end:


#Recip(HG): the reciprocal 2F1 of HG
Recip:=proc(HG) :
if not (nops(HG)=3 and nops(HG[1])=2 and nops(HG[2])=1) then
 print(`Not a 2F1`):
 RETURN(FAIL):
fi:

[[HG[1][1],HG[1][1]-HG[2][1]+1],[HG[1][1]-HG[1][2]+1],1/HG[3]]:
end:


#PfaffT2(HG): Applies (2.2.6) in AAR
PfaffT2:=proc(HG):
if not (nops(HG)=3 and nops(HG[1])=2 and nops(HG[2])=1) then
 print(`Not a 2F1`):
 RETURN(FAIL):
fi:
[[HG[1][1], HG[2][1]-HG[1][2]],HG[2],HG[3]/(HG[3]-1)]:end:

#Euler(HG): Applies (2.2.7) in AAR
Euler:=proc(HG):
if not (nops(HG)=3 and nops(HG[1])=2 and nops(HG[2])=1) then
 print(`Not a 2F1`):
 RETURN(FAIL):
fi:
[[HG[2][1]-HG[1][1], HG[2][1]-HG[1][2]],HG[2],HG[3]]:end:

#Quad1(HG): Applies (3.1.3) in AAR
Quad1:=proc(HG):
if not (nops(HG)=3 and nops(HG[1])=2 and nops(HG[2])=1) then
 print(`Not a 2F1`):
 RETURN(FAIL):
fi:

if 2*HG[2][1]-HG[1][1]-HG[1][2]-1<>0 then
print(`Not applicable`):

RETURN(FAIL):
fi:

[[HG[1][1]/2,HG[1][2]/2],HG[2],4*HG[3]*(1-HG[3])]:end:

#Quad2(HG): Applies (3.1.7) in AAR
Quad2:=proc(HG):
if not (nops(HG)=3 and nops(HG[1])=2 and nops(HG[2])=1) then
 print(`Not a 2F1`):
 RETURN(FAIL):
fi:

if 2*HG[1][1]-HG[2][1]<>0 then
print(`Not applicable`):
RETURN(FAIL):
fi:

[[HG[1][2]/2,(HG[1][2]+1)/2],HG[1][1]+1/2, (HG[3]/(2-HG[3]))^2]:end:



#Khaverim1(tri1): given a triple tri1=[A,B,C], finds
#all the immediately
#implied triples by Neg. Symmm. Recipr. Pfaff2, Euler and Pfaff1
Khaverim1:=proc(tri1)
local A,B,C:
 A:=tri1[1]: B:=tri1[2]: C:=tri1[3]:

{[-A,-B,-C],[B,A,C],[A,-C+A,-B+A],[A,C-B,C],[C-A,C-B,C],
[A,B,B+1+A-C]}:

end:

#Khaverim1S(tri1S): the set of friends of all the members of tri1S
Khaverim1S:=proc(tri1S)
local i:
{ seq(op(Khaverim1(tri1S[i]) ),i=1..nops(tri1S))}:
end:


#IKhaverimS1(tri1,L): all the iterated friends up to L times
IKhaverimS1:=proc(tri1,L)
local gu,i:
gu:={tri1}:

for i from 1 to L do
gu:=gu union Khaverim1S(gu):
od:
gu:
end:


#IKhaverimS(tri1S,L): the set of friends (up to
#L degree os separation) of all the members of tri1S
IKhaverimS:=proc(tri1S,L)
local i:
{ seq(op(IKhaverimS1(tri1S[i],L) ),i=1..nops(tri1S))}:
end:

#CandSet(K,L): a;; the triples of integers between
#-K and K that are inequivalent to each other up to
#L degrees of separation
CandSet:=proc(K,L)
local gu,i1,i2,i3,lu:
gu:={seq(seq(seq([-i1,i2,i3],i1=1..K),i2=-K..K),i3=-K..K)}:

lu:={}:


while gu<>{} do
lu:=lu union {gu[1]}:
gu:=gu minus IKhaverimS1(gu[1],L):
od:

lu:

end:



#FindHGr1(tri1,n,R): FindHGr([[-a*n,b*n+b1],[c*n+c1],z] ,n,{b1,c1,z},R):
FindHGr1:=proc(tri1,n,R) local a,b,c,b1,c1,z: a:=tri1[1]: b:=tri1[2]: c:=tri1[3]:
FindHGr([[a*n,b*n+b1],[c*n+c1],z] ,n,{b1,c1,z},R):
end:




#DoronD(F,k,n,N): Does Doron with F given in hypergeometric notation [Numer,Den,z]
DoronD:=proc(F,k,n,N)
Doron(HGtoTerm(F,k),k,n,N):
end:


#Implications(HG): the set of implications of HG by Recip, Euler, PfaffT1,and PfaffT2
Implications:=proc(HG):
{[[HG[1][2],HG[1][1]],HG[2],HG[3]],Recip(HG),Euler(HG),PfaffT1(HG),PfaffT2(HG)}:
end:


#ImplicationsL1(HG,L): all the implications of an exact evaluation HG with exactly L 
ImplicationsL1:=proc(HG,L) local gu,i1:
if L=1 then
 RETURN(Implications(HG)):
fi:

gu:=ImplicationsL(HG,L-1):

{seq(op(Implications(gu[i1])),i1=1..nops(gu))}:
end:

#ImplicationsL(HG,L):the implications of an exact evaluation HG up
ImplicationsL:=proc(HG,L) local i1:
{seq(op(ImplicationsL1(HG,i1)),i1=1..L)}:
end:

#Find2F1R(K,n,R,L): 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
#and such that none follows from others by one of the rules
#in Implications, up to L iterations
#try Find2F1R(2,n,R)
Find2F1R:=proc(K,n,R,L)
local T,gu,mu,b1,c1,z,mu1,a,b,c,gu1,i,ku,i1,ku1,i2,i3,muam,gu2:
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) and
        not member(mu1[i][1], {seq(op(ImplicationsL(gu1[i1],L)),i1=1..nops(gu1))}) 
        and NotRedundant(mu1[i][1]) and not IsKnown(mu1[i][1]) 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: 
 
ku1:=[]:

for i1 from 1 to K do
 for i2 from -K to K do
  for i3 from -K to K do
   muam:=[-i1,i2,i3]:
    if member(muam, ku) then
     ku1:=[op(ku1),muam]:
    fi:
  od:
 od:
od:

gu2:=[]:
for i1 from 1 to nops(ku1) do
 gu2:=[op(gu2),op(T[op(ku1[i1])])]:
od:


ku1,op(T),gu2:
end:

NotRedundant:=proc(HG)
local i,j:
for i from 1 to nops(HG[1])  do

  if type(HG[1][i],integer) and  HG[1][i]>=0 then
    RETURN(false):
  fi:

 for j from 1 to nops(HG[2])  do
  if type(HG[1][i]-HG[2][j],integer) and  HG[1][i]-HG[2][j]>=0 then
    RETURN(false):
  fi:

 od:
od:
true:
end:



IsKnown1:=proc(HG):
if HG[3]=1 then
 RETURN(true):
elif HG[3]=-1 and HG[1][1]-HG[1][2]+1=HG[2][1] then
  RETURN(true):
elif HG[3]=-1 and HG[1][2]-HG[1][1]+1=HG[2][1] then
  RETURN(true):
elif HG[3]=1/2 and HG[1][1]+HG[1][2]=1 then
   RETURN(true):
elif HG[3]=1/2 and HG[1][1]+HG[1][2]+1=2*HG[2][1] then
   RETURN(true):
else
  RETURN(false):
 fi:
end:

IsKnown:=proc(HG)
local lu,i:
lu:={HG} union Implications(HG):
for i from 1 to nops(lu) do
 if IsKnown1(lu[i]) then
  RETURN(true):
 fi:
od:
false:
end:


Yafe:=proc(i,HG) local lu,mu:
lu:=HG[1]: mu:=HG[2]:
printf("Theorem %a : F(%a , %a , %a , %a )=%a .\n", i, lu[1][1],lu[1][2],lu[2][1],lu[3] , mu);

printf("\n"):
printf("Proof : %a . QED.",HG[3]):
printf("\n"):
printf("\n"):
end:


#Sefer(S): writes a book using the info in S
Sefer:=proc(S) local i:
print(`This book contains`, nops(S), `statements and proofs of exact 2F1 hypergeometric evaluations`):
print(``):
print(`Let R(a,k):=a(a+1)...(a+k-1)`):
for i from 1 to nops(S) do
 Yafe(i,S[i]):
od:
end:

Book2F1:=proc(K,n,R,L)
 Sefer(Find2F1R(K,n,R,L)[3]):
end:

BOOK:=proc(K,K1,n,R,L)
 Sefer(Find2F1Ra(K,K1,n,R,L)[3]):
end:



#Find2F1Ra(K,K1,n,R,L): 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<=K1, -K<= b <= K, -K <= c <= K c<>{a,b}
#and not consequences of previous ones
#R is the symbol for rising-factorial
#and such that none follows from others by one of the rules
#in Implications, up to L iterations
#try Find2F1Ra(2,n,R)
Find2F1Ra:=proc(K,K1,n,R,L)
local T,gu,mu,b1,c1,z,mu1,a,b,c,gu1,i,ku,i1,ku1,i2,i3,muam,gu2,b2,c2:
ku:={}:
gu:={}:
gu1:={}:
for a from 1 to K1 do
for b2 from 0 to K do
 for c2 from 0 to K do

   b:=b2: c:=c2:
  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) and
        not member(mu1[i][1], {seq(op(ImplicationsL(gu1[i1],L)),i1=1..nops(gu1))}) 
        and NotRedundant(mu1[i][1]) and not IsKnown(mu1[i][1]) 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:



   b:=b2: c:=-c2:
  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) and
        not member(mu1[i][1], {seq(op(ImplicationsL(gu1[i1],L)),i1=1..nops(gu1))}) 
        and NotRedundant(mu1[i][1]) and not IsKnown(mu1[i][1]) 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:
  
  
   b:=-b2: c:=c2:
  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) and
        not member(mu1[i][1], {seq(op(ImplicationsL(gu1[i1],L)),i1=1..nops(gu1))}) 
        and NotRedundant(mu1[i][1]) and not IsKnown(mu1[i][1]) 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:
  
   b:=-b2: c:=-c2:
  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) and
        not member(mu1[i][1], {seq(op(ImplicationsL(gu1[i1],L)),i1=1..nops(gu1))}) 
        and NotRedundant(mu1[i][1]) and not IsKnown(mu1[i][1]) 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: 

ku1:=[]:

for i1 from 1 to K do
 for i2 from -K to K do
  for i3 from -K to K do
   muam:=[-i1,i2,i3]:
    if member(muam, ku) then
     ku1:=[op(ku1),muam]:
    fi:
  od:
 od:
od:

gu2:=[]:
for i1 from 1 to nops(ku1) do
 gu2:=[op(gu2),op(T[op(ku1[i1])])]:
od:

ku1,op(T),gu2:
 

end:


#IsMultiple(HG,n): Is it the case n=r*n for some r
IsMultiple:=proc(HG,n)
local gu,HG1,k,N:
gu:=gcd(gcd(coeff(HG[1][1],n,1), coeff(HG[1][2],n,1)), coeff(HG[2][1],n,1)):

if gu=1 then
  RETURN(false):
fi:

HG1:=subs(n=n/gu,HG):

print(Doron(HGtoTerm(HG1,k),k,n,N)[1]):
print(degree(Doron(HGtoTerm(HG1,k)[1],k,n,N)[1],N)):

if degree(Doron(HGtoTerm(HG1,k),k,n,N)[1],N)=1 then
 RETURN(true):
else
 RETURN(false):
fi:

end:


#HypToInt1(Hyp,y): converting a 2F1 to integral form
HypToInt1:=proc(Hyp,y)
local a,b,c,z:
a:=Hyp[1][1]: b:=Hyp[1][2]: c:=Hyp[2][1]: z:=Hyp[3]:
y^(b-1)*(1-y)^(c-b-1)*(1-y*z)^(-a)*(c-1)!/(b-1)!/(c-b-1)!:
end:

#HypToInt(Hyp,Rhs,y): converting a 2F1 Hyp to integral form
#with Rhs given in terms of R(a,k) denoting rising factiroal
#For example, try:
#HypToInt([[-n,a],[b],1],(c-1)!*(c+n-b-1)!/(c+n-1)!/(c-b-1)!,y);
#
HypToInt:=proc(Hyp,Rhs,y)
local a,b,c,z:
a:=Hyp[1][1]: b:=Hyp[1][2]: c:=Hyp[2][1]: z:=Hyp[3]:
y^(b-1)*(1-y)^(c-b-1)*(1-y*z)^(-a)*(c-1)!/(b-1)!/(c-b-1)!/
subs(R=Rf1,Rhs):
end:
 
####From EKHAD8



goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:

duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1:
KAMA:=10:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:

pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:

 
AZd:= proc(A,y,n,N) local ORDER,gu,KAMA:
KAMA:=8:
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:

###End portion from EKHAD8

 

rf1:=proc(a,k):
if a=0 then 
RETURN(0):
else
(a+k-1)!/(a-1)!
fi:
end:

#FindHGc(HG,n,Pars,R,y): Discrete-Continuous WZ-pairs implied by FindHGr
FindHGc:=proc(HG,n,Pars,y) local mu,gu,i,F,ku,N,coe:
mu:=FindHGr(HG,n,Pars,rf1) :

gu:={}:

for  i from 1 to nops(mu) do
 F:=HypToInt(mu[i][1],mu[i][2],y):
 ku:=AZd(F,y,n,N): 
 coe:=normal(ku[1]/(N-1)):
 if not type(coe,integer) then
  RETURN(FAIL):
 fi:
gu:= gu union { [normal(F/coe),normal(ku[2]*F) ] }:
od:
gu:
end:


Find2F1wz:=proc(K,K1,n,L,y)
local gu,mu,i,ku,N,coe,F,R:
gu:={}:
mu:=Find2F1Ra(K,K1,n,R,L)[3]:
mu:=subs(R=rf1,mu):

for  i from 1 to nops(mu) do
 F:=HypToInt(mu[i][1],mu[i][2],y):
 ku:=AZd(F,y,n,N): 
 coe:=normal(ku[1]/(N-1)):
gu:= gu union { [F,normal(ku[2]*F/coe) ] }:
od:
gu:

end: