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

read EKHAD:
print(` Version of June 13, 2002`):
print(`First version, ca. Sept. 12, 2001`):

print(`This is AperyWZ, a Maple package that tries to find`):
print(`Apery-type irrationality proofs.`):
print(` using WZ pairs, as initially described in`):
print(` Closed Form (pun intended!), Conemporary Mathematics v. 143`):
print(` edited by M. Knopp and Mark Sheingorn, AMS, 1993, 579-608 `):
print(` available from http://www.math.rutgers.edu/~zeilberg/ `):
print(``):
print(`It is one of the packages acompanying Zeilberger's article`):
print(`COMPUTERIZED DECONSTRUCTION`):
print(` It requires EKHAD in the same directory`):

print(`The most current version `):
print(`is available at http://www.math.rutgers.edu/~zeilberg.`):
print(`Please report all bugs to: zeilberg@math.rutgers.edu .`):
lprint(``):
print(`Written by Doron Zeilberger`):
print(``):
print(`For general help type: Ezra();`):
 
Ezra:=proc()
if args=NULL then
print(`This is AperyWZ, a Maple packages that tries to find`):
print(`Apery-type irrationality proofs.`):
print(` using WZ pairs`):
print(`The most current version of the program and article `):
print(`are available at http://www.math.rutgers.edu/~zeilberg.`):
print(`Please report all bugs to: zeilberg@math.rutgers.edu .`):
lprint(``):
print(`Written by Doron Zeilberger`):
lprint(``):
print(` Contains Procedures: `):
print(``):
print(`Acc, an`):
print(` beta , bn, CheckAcc, CheckWZ, cn, cnk, cnkS, Constant1 `):
print(` delta, deltaWZb, Natcn, WZchu, WZdixon, WZkummer,`):
print(` WZpair, yakhas, zeilWZP, zeilWZPv `):
print(` `):
print(`To get a list of built in WZ-pairs, type: Ezra(WZdata):`):
print(``):
print(` For specific help type "Ezra(procedure_name)" `):
fi:

if nops([args])=1 and op(1,[args])=Acc then
print(`Acc(F,G,k,n): the acceleration formula implied by the WZ`):
print(`pair (F,G) in the variables k,n`):
print(`For example try Acc(WZpair((-1)^*k*(n+k)!/k!^2/(n-k)!,k,n),k,n);`):
fi:

if nops([args])=1 and op(1,[args])=an then
print(`an(F,G,bnk,n,k,n0): Given a WZ pair: F,G  and a weighted`):
print(`averaging function, bnk (an expression in n and k)`):
print(`finds the sum(cnk(n0,k)bnk(n0,k), k=0..n0): For example, `):
print(`try an(WZlog2,binomial(n,k)*binomial(n+k,k),n,k,5);`):
print(`try an(WZzeta2,binomial(n,k)^2*binomial(n+k,k),n,k,5);`):
print(`try an(WZzeta3,binomial(n,k)^2*binomial(n+k,k)^2,n,k,5);`):
fi:

if nops([args])=1 and op(1,[args])=beta then
print(`beta(ope1,N,n): Given a recurrence operator, ope1, in the`):
print(`shift operator N, and the variable n, finds the beta such`):
print(`that the sequences satisfying the recurrence are O((beta)^n)`):
fi:

if nops([args])=1 and op(1,[args])=bn then
print(` bn(bnk,n,k,n0): Given an `):
print(` averaging function, bnk (an expression in n and k) `):
print(` finds the sum(bnk(n0,k), k=0..n0): For example, `):
print(` try bn(binomial(n,k)^2*binomial(n+k,k),n,k,5); `):
fi:

if nops([args])=1 and op(1,[args])=CheckAcc then
print(`CheckAcc(R1,R2,n,N): finds the difference between the sum of`):
print(` R2(n) from n=1 to N and R1(n) from n=0 to N `):
fi:

if nops([args])=1 and op(1,[args])=CheckWZ then
print(`CheckWZ(F,G,k,n): checks whether (F,G) is a WZ pair`):
print(`in the variables k,n`):
fi:

if nops([args])=1 and op(1,[args])=cn then
print(` cn(F,G,bnk,n,k,n0): Given a WZ-pair, and  `):
print( ` averaging function, bnk (all expressions in n and k) `):
print(` and an integer  n0, finds the rational number `):
print(` an(F,G,bnk,n,k,n0):/bn(bnk,n,k,n0): `):
print(` finds the sum(bnk(n0,k), k=0..n0): For example,  `):
print(` try cn(WZlog2, binomial(n,k)^2*binomial(n+k,k),n,k,5);`):
fi:

if nops([args])=1 and op(1,[args])=cnk then
print(`cnk(F,G,n,k,n0,k0): Given a WZ pair F,G, in the discrete`):
print(`	 variables n,k, and specific integers n0,k0`):
print(`Finds the potential function c(n,k) (such that F= DELTA_k c`):
print(` and G= DELTA_n c, at n=n0,k=k0`):
print(` For example try cnk(WZzeta2,n,k,5,3)`):
fi:

if nops([args])=1 and op(1,[args])=cnkS then
print(`cnkS(F,G,n,k,m): Given a WZ pair F,G, in the discrete variables`):
print(`m, and a (discrete) variable m returns the sumbolic expression`):
print(`for the potential function c(n,k) (such that F= DELTA_k c`):
print(` and G= DELTA_n c, using the summation variable m`):
print(`For example, try: cnkS(WZzeta2,n,k,m); `):
fi:


if nops([args])=1 and op(1,[args])=Constant1 then
print(` Constant1(R1,R2,n,N): Given an acclerating pair `):
print(` finds the constant up to error 1/10^N in the  `):
print( `series term `):
fi:

if nops([args])=1 and op(1,[args])=delta then
print(`delta(sid,K): Given a sequence of rational numbers, sid`):
print(` and a real number K `):
print(` finds the smallest delta in the list s.t. `):
print(` |K-an/bn|<=C/bn^(1+delta) `):
print(` where K is its estiamted limit, the last term `):
print(` For example, try: with(combiant): `):
print(` delta([seq(fibonacci(i+1)/fibonacci(i),i=1..100)]): `):
fi:

if nops([args])=1 and op(1,[args])=deltaWZb  then
print(`deltaWZb(F,G,b,n,k,K,N): Given a WZ pair F,G, and a weighting`):
print(`function b, all in terms of the discrete variables n and k`):
print(`and a real constant K, and an integer N finds, empirically`):
print(` the "irattionality factor" delta, for the first N terms `):
print(`of the sequence cn(F,G,bnk,n,k,n0), of approximants to K`):
print(`according to Apery's method`):
print(` For example try: `):
print(`deltaWZb(WZlog2,binomial(n,k)*binomial(n+k,k),n,k, `):
print( `evalf(log(2)),30);` ):
print(`deltaWZb(WZzeta2,binomial(n,k)^2*binomial(n+k,k), `):
print(` n,k,evalf(Pi^2/6),30);`):
print(` deltaWZb(WZzeta3,binomial(n,k)^2*binomial(n+k,k)^2,n,k, `):
print(` evalf(Zeta(3)),30); `):
fi:

if nops([args])=1 and op(1,[args])=WZchu then
print(`WZchu(n,k,a,b,c,bitui): `):
print(` The generalizations of the WZ pair for Zeta(2)`):
print(`obtained by shadowing the WZ-pair implied bu Vandermonde-Chu`):
print(`(alias Gauss's 2F1(1)), bitui is an expressions in n`):
print(` a,b,c, are numbers, parameters, or expressions in n`):
print(`For example, try: WZchu(n,k,0,0,0,n);`):
fi:

if nops([args])=1 and op(1,[args])=WZdixon then
print(`WZdixon(n,k,a,c,bitui): The generalizations of the WZ pair`):
print(` for Zeta(3)`):
print(`obtained by shadowing the WZ-pair implied by Dixon's identity`):
print(`with subsitituting to k bitui (starting with G)`):
print(`For example, try: WZdixon(n,k,0,0,k)`):
fi:

if nops([args])=1 and op(1,[args])=Natcn then
print(` Natcn(F,G,n,k,n0): Given a WZ pair: F,G   `):
print(` ( expressions in n and k) `):
print(` finds the sum(cnk(n0,k)*denom(cnk(n0,k0)), k=0..n0): `):
print(` deivided by sum(denom(cnk(n0,k0)), k=0..n0):  `):
print(` For example, try Natcn(WZzeta2,n,k,5); `):
fi:

if nops([args])=1 and op(1,[args])=WZdata then
print(` The built-in WZ pairs are: WZlog1x, WZlog2,WZkummer0, `):
print(` WZzeta2, WZzeta3`):
fi:


if nops([args])=1 and op(1,[args])=WZkummer then
print(`WZkummer(n,k,b,bitui): A WZ pair obtained from Kummer's 2F1(-1)`):
print(`identity, obtained by subsitituting to n bitui and  b is a number`):
print(`or paremeter, For example, try: WZkummer(n,k,0,n); `):
fi:

if nops([args])=1 and op(1,[args])=WZpair then
print(`WZpair(F,k,n): Gives the WZ pair induced by the closed form`):
print(`function F (in variables n and k), For example, try `):
print(`WZpair(binomial(n,k),k,n):`):
fi:

if nops([args])=1 and op(1,[args])=yakhas then
print(`yakhas(hg,n): the asymptoic ratio of hg(n+1)/hg(n)`):
fi:

if nops([args])=1 and op(1,[args])=zeilWZP then
print(`Note: this is a terse version, for the terse version`):
print(`Use zeilWZPv `):
print(`zeilWZP(bnk,F,G,k,n,N), given an averaging function bnk`):
print(`and a WZ pair, (F,G) (all expressions in n,k) and a shift`):
print(`operator N, finds the recurrences satisfied by the`):
print(`numerator and denominators in the Apery-style approximation`):
print(`and the certificates, in the form: ope1,cert1,ope2,cert2`):
print(``):
print(`For example, try: `):
print(` zeilWZP(binomial(n,k)*binomial(n+k,k),WZlog2,k,n,N) `):
print(` zeilWZP(binomial(n,k)^2*binomial(n+k,k),WZzeta2,k,n,N) `):
print(` zeilWZP(binomial(n,k)^2*binomial(n+k,k)^2,Wzeta3,k,n,N) `):

fi:

if nops([args])=1 and op(1,[args])=zeilWZPv then
print(`Note: this is a verbose version, for the terse version`):
print(`Use zeilWZP `):
print(`zeilWZPv(bnk,F,G,k,n,N), given an averaging function bnk`):
print(`and a WZ pair, (F,G) (all expressions in n,k) and a shift`):
print(`operator N, finds the recurrences satisfied by the`):
print(`numerator and denominators in the Apery-style approximation`):
print(``):
print(`For example, try: `):
print(` zeilWZPv(binomial(n,k)*binomial(n+k,k),WZlog2,k,n,N) `):
print(` zeilWZPv(binomial(n,k)^2*binomial(n+k,k),WZzeta2,k,n,N) `):
print(` zeilWZPv(binomial(n,k)^2*binomial(n+k,k)^2,Wzeta3,k,n,N) `):

fi:


end:

####A data-base of WZ pairs
WZlog1x:=(1+x)*(-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1),
x*(-1)^k*n!*k!/(n+k+1)!/x^k/(1+x)^(n+1):

WZlog2:=(2)*(-1)^k*n!*k!/(n+k+1)!/2^(n+1),
(-1)^k*n!*k!/(n+k+1)!/2^(n+1):

WZzeta2:=(-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!,
(-1)^(n+k)*k!^2*(n-k)!/(n+k+1)!*2/(n+1):


WZzeta3:=(-1)^k*k!^2*(n-k-1)!/(n+k+1)!/2/(k+1),
(-1)^k*k!^2*(n-k)!/(n+k+1)!/(n+1)^2:

WZkummer0:=
(-1)^k*k!*(k-b)!*(2*n+b)!*(n+b)!*n!/b!^2/(2*n+b+k+1)/(2*n+b+k)!/(2*n+k+1)/
(2*n+k)!/((-1)^n) ,  
1/2*(17*n+b^2+10*n^2+6*b+7+7*n*b+2*k*b+6*n*k+5*k+k^2)*(-1)^k*k!*
(k-b)!*(2*n+b)!*(n+b)!*n!/b!^2/(2*n+b+k+2)/
(2*n+b+k+1)/(2*n+b+k)!/(2*n+k+2)/(2*n+k+1)/(2*n+k)!/((-1)^n):

###

#beta(ope1,N,n): Given a recurrence operator, ope1, in the
#shift operator N, and the variable n, finds the beta such
#that the sequences satisfying the recurrence are O((beta)^n)
beta:=proc(ope1,N,n)
local gu,deg,gu1:
gu:=expand(ope1):
deg:=degree(gu,n):
gu1:=coeff(gu,n,deg):
solve(gu1,N):
end:


#WZpair(F,k,n): Gives the WZ pair induced by the closed form
#function F (in variables n and k), For example, try 
#WZpair(binomial(n,k),n,k):
WZpair:=proc(F,k,n)
local F1,G1,ope,N,r,mu,lu,ku:

ope:=expand(zeil(F,k,n,N)[1]):

if degree(ope,N)<>1 then
 print(`Operator not first order`):
 RETURN(0):
fi:

mu:=rsolve({coeff(ope,N,0)*r(n)+coeff(ope,N,1)*r(n+1),r(0)=1},r):
F1:=simplify(F/mu):
lu:=zeil(F1,k,n,N):
ku:=normal(lu[1]/(N-1)):

if diff(ku,N)<>0 then
ERROR(`Something is wrong`):
fi:
G1:=lu[2]*F1/ku:
F1:=convert(F1,factorial):
G1:=convert(G1,factorial):
F1:=normal(expand(F1)):
G1:=normal(expand(G1)):
F1,G1:
#simplify(F1),simplify(G1):

end:

#CheckWZ(F,G,k,n): checks whether (F,G) is a WZ pair
#in the variables k,n
CheckWZ:=proc(F,G,k,n):
normal(simplify((subs(n=n+1,F)-F-(subs(k=k+1,G)-G))/F)):
end:


#CheckWZS(F,G,k,n,k0,n0): checks whether (F,G) is a WZ pair
#in the variables k,n at the point (n0,k0)
CheckWZS:=proc(F,G,k,n,k0,n0):
subs({k=k0,n=n0+1},F)-
subs({n=n0,k=k0},F)-
subs({n=n0,k=k0+1},G)+
subs({n=n0,k=k0},G):
end:

#Acc(F,G,k,n): the acceleration formula implied by the WZ
#pair (F,G) in the variables k,n
Acc:=proc(F,G,k,n):

if CheckWZ(F,G,k,n)<>0 then
 ERROR(`Input is not a WZ-pair`):
fi:

normal(simplify(subs({k=0},G))),
normal(simplify(subs(k=n-1,F)+subs({k=n-1,n=n-1},simplify(G)))):
end:

#CheckAcc(R1,R2,n,N): finds the difference between the sum of
#R2(n) from n=1 to N and R1(n) from n=0 to N
CheckAcc:=proc(R1,R2,n,N)
local i,mu1,mu2:

mu1:=evalf(subs(n=0,R1)):
mu2:=0:

for i from 1 to N do
mu1:=mu1+evalf(subs(n=i,R1)):
mu2:=mu2+evalf(subs(n=i,R2)):
od:

mu1,mu2,mu1-mu2:

end:

#cnkS(F,G,n,k,m): Given a WZ pair F,G, in the discrete variables
#n,k, and a (discrete) variable m returns the sumbolic expression
#for the potential function c(n,k) (such that F= DELTA_k c
# and G= DELTA_n c, using the summation variable m
cnkS:=proc(F,G,n,k,m):
Sum(normal(simplify(subs({k=0,n=m-1},G))),m=1..n)+
Sum(normal(simplify(subs({k=m-1},F))),m=1..k):
end:


#cnkOld(F,G,n,k,n0,k0): Given a WZ pair F,G, in the discrete variables
#n,k, and specific integers n0,k0
#Finds the potential function c(n,k) (such that F= DELTA_k c
# and G= DELTA_n c, at n=n0,k=k0
# For example try cnk(WZzeta2,n,k,5,3)
cnkOld:=proc(F,G,n,k,n0,k0) local m:
value(Sum(normal(simplify(subs({k=0,n=m-1},G))),m=1..n0))+
value(Sum(normal(simplify(subs({k=m-1,n=n0},F))),m=1..k0)):
end:

cnkOld:=proc(F,G,n,k,n0,k0) local m,lu:
lu:=cnkS(F,G,n,k,m):
value(convert(subs({n=n0,k=k0},lu),factorial)):
end:


#an(F,G,bnk,n,k,n0): Given a WZ pair: F,G  and a weighted
#averaging function, bnk (an expression in n and k)
#finds the sum(cnk(n0,k)bnk(n0,k), k=0..n0): For example, 
#try an(WZzeta2,binomial(n,k)^2*binomial(n+k,k),n,k,5);
an:=proc(F,G,bnk,n,k,n0)
local gu,k0:
gu:=0:
for k0 from 0 to n0 do
gu:=gu+eval(subs({n=n0,k=k0},bnk)*cnk(F,G,n,k,n0,k0)):
od:
gu:
end:


#bn(bnk,n,k,n0): Given an
#averaging function, bnk (an expression in n and k)
#finds the sum(bnk(n0,k), k=0..n0): For example, 
#try bn(binomial(n,k)^2*binomial(n+k,k),n,k,5);
bn:=proc(bnk,n,k,n0)
local gu,k0:
gu:=0:
for k0 from 0 to n0 do
gu:=gu+eval(subs({n=n0,k=k0},bnk)):
od:
gu:
end:


#cn(F,G,bnk,n,k,n0): Given a WZ-pair, and 
#averaging function, bnk (all expressions in n and k)
#and an integer  n0, finds the rational number
#an(F,G,bnk,n,k,n0):/bn(bnk,n,k,n0):
#finds the sum(bnk(n0,k), k=0..n0): For example, 
#try cn(WZlog2, binomial(n,k)^2*binomial(n+k,k),n,k,5);
cn:=proc(F,G,bnk,n,k,n0):an(F,G,bnk,n,k,n0)/bn(bnk,n,k,n0):end:


#delta(sid,K): Given a sequence of rational numbers, sid
#and a real number K
#finds the smallest delta in the list s.t. |K-an/bn|<=C/bn^(1+delta)
#where K is its estiamted limit, the last term
#For example, try: with(combiant):
#delta([seq(fibonacci(i+1)/fibonacci(i),i=1..100)]):
delta:=proc(sid,K)
local n,mu,i:

n:=nops(sid):

mu:=[]:
for i from 1 to n do
mu:=[op(mu),evalf(-log(abs(K-sid[i]))/log(denom(sid[i])))-1]:
od:
min(op(mu)):
end:


#deltaWZb(F,G,b,n,k,K,N): Given a WZ pair F,G, and a weighting
#function b, all in terms of the discrete variables n and k
#and a real constant K, and an integer N finds, empirically
#the "irattionality factor" delta, for the first N terms
#of the sequence cn(F,G,bnk,n,k,n0), of approximants to K
#according to Apery's method
#For example try:
#deltaWZb(WZlog2,binomial(n,k)*binomial(n+k,k),n,k,evalf(log(2)),30);
#deltaWZb(WZzeta2,binomial(n,k)^2*binomial(n+k,k),n,k,evalf(Pi^2/6),30);
#deltaWZb(WZzeta3,binomial(n,k)^2*binomial(n+k,k)^2,n,k,
#evalf(Zeta(3)),30);
deltaWZb:=proc(F,G,b,n,k,K,N)
local lu,n0:
lu:=[seq(cn(F,G,b,n,k,n0),n0=1..N)]:
delta(lu,K):
end:


#Natcn(F,G,n,k,n0): Given a WZ pair: F,G  
#( expressions in n and k)
#finds the sum(cnk(n0,k)*denom(cnk(n0,k0)), k=0..n0): 
#deivided by sum(denom(cnk(n0,k0)), k=0..n0): 
#For example, try Natcn(WZzeta2,n,k,5);
Natcn:=proc(F,G,n,k,n0)
local gu1,gu2,k0,mu:
gu1:=0: gu2:=0:
for k0 from 0 to n0 do
 mu:=cnk(F,G,n,k,n0,k0):
 gu1:=gu1+mu*denom(mu):
 gu2:=gu2+denom(mu):
od:
gu1/gu2:
end:

 
#zeilWZPv(bnk,F,G,k,n,N), given an averaging function bnk
#and a WZ pair, (F,G) (all expressions in n,k) and a shift
#operator N, finds the recurrences satisfied by the
#numerator and denominators in the Apery-style approximation
#
zeilWZPv:=proc(bnk,F,G,k,n,N)
local i,j,gu,lu,R,cert,ku:

ku:=simplify(F/G):

if not type(ku,algfun(rational)) then
ERROR(` the second and third arguments are not WZ `):
fi:

gu:=simplify(subs(n=n+1,F)/F)-1-simplify(subs(k=k+1,G)/F)+simplify(G/F):

gu:=normal(gu):

if gu<>0  then
ERROR(` the second and third arguments are not WZ `):
fi:


gu:=zeil(bnk,k,n,N):

R:=gu[1]:
cert:=gu[2]:

gu:=0:

for i from 0 to degree(R,N) do
lu:=0:
for j from 0 to i-1 do
lu:=lu+normal(simplify(subs(n=n+j,G)/F)):
lu:=normal(lu):
od:

gu:=gu+coeff(R,N,i)*simplify1(bnk,n,i)*lu:
od:
gu:=gu-subs(k=k+1,cert)*simplify1(bnk,k,1):

gu:=normal(gu):
gu:=gu*bnk*F:


#print(`H(n,k) is`,gu):
gu:=zeil(gu,k,n,N):
print(`Theorem: Let b(`,n,k,`) be defined by`):
print(bnk):
print(`and let c(`,n,k,`) be the potential function of the WZ pair (F,G)`):
print(`where F and G are `):
print(F,G):
print(`Let a(`,n,`) be the sum of F(`,n,k,`)*c(`,n,k,`) `):
print(`this sum is annihilated by the operator`):
print(gu[1], `times` , R):
lprint(``):
print(`Proof by Shalosh B. Ekhad`):
print(`We first find the linear recurrence operator annihilating`):
print(`the sum of b(`,n,k,`), which is `):
print(R):
print(`For a proof do zeil(b,k,n) in EKHAD.`):
print(`The certificate turned out to be`):
print(cert):
print(`Applying this operator to c(n,k)*b(n,k) and subtrating`):
print(`c(n,k) times the operator applied on b(n,k), as was done`):
print(`in Doron Zeilberger's article: Closed Form (pun intended!)`):
print(`Contemporary Mathematics 143, pp 579-607, Theorem 9 (p. 600 ff) `):
print(`with slight notational changes to accomodate EKHAD's current`):
print(`convention of using forward differnce operators.`):
print(`We get that the operator S(N,n) promised by Theorem 9, is`):
print(gu[1]):
print(`The certificate certifying S(N,n), i.e. the one defining`):
print(`D(n,k), was found by EKHAD to be`):
print(gu[2]):


R,cert,gu:

end:
 
 
#zeilWZP(bnk,F,G,k,n,N), given an averaging function bnk
#and a WZ pair, (F,G) (all expressions in n,k) and a shift
#operator N, finds the recurrences satisfied by the
#numerator and denominators in the Apery-style approximation
#and the certificates, terse style
zeilWZP:=proc(bnk,F,G,k,n,N)
local i,j,gu,lu,R,cert,ku:

ku:=simplify(F/G):

if not type(ku,algfun(rational)) then
ERROR(` the second and third arguments are not WZ `):
fi:

gu:=simplify(subs(n=n+1,F)/F)-1-simplify(subs(k=k+1,G)/F)+simplify(G/F):

gu:=normal(gu):

if gu<>0  then
ERROR(` the second and third arguments are not WZ `):
fi:


gu:=zeil(bnk,k,n,N):

R:=gu[1]:
cert:=gu[2]:

gu:=0:

for i from 0 to degree(R,N) do
lu:=0:
for j from 0 to i-1 do
lu:=lu+normal(simplify(subs(n=n+j,G)/F)):
lu:=normal(lu):
od:

gu:=gu+coeff(R,N,i)*simplify1(bnk,n,i)*lu:
od:
gu:=gu-subs(k=k+1,cert)*simplify1(bnk,k,1):

gu:=normal(gu):
gu:=gu*bnk*F:


gu:=zeil(gu,k,n,N):

R,cert,gu:

end:


#GWZzeta2(n,k,a,b,c): The generalizations of the WZ pair for Zeta(2)
#obtained by shadowing the WZ-pair implied bu Vandermonde-Chu
#(alias Gauss's 2F1(1)), a,b,c are the parameters
GWZzeta2:=proc(n,k,a,b,c)
local F:
F:=(-1)^k*k!*(k+c)!*(n-k-1+a)!/(n+k+1+b)!:
WZpair(F,k,n):
end:

#WZchu(n,k,a,b,c,bitui): The generalizations of the WZ pair for Zeta(2)
#obtained by shadowing the WZ-pair implied by Vandermonde-Chu
#with subsitituting to n bitui
#(alias Gauss's 2F1(1)), a,b,c are the parameters
WZchu:=proc(n,k,a,b,c,bitui)
local F:
F:=subs(n=bitui,(-1)^k*k!*(k+c)!*(n-k-1+a)!/(n+k+1+b)!):

WZpair(F,k,n):
end:


cnk:=proc(F,G,n,k,n0,k0) local m,gu:
gu:=0:

for m from 0 to n0-1 do
gu:=gu+subs({k=0,n=m},G):
od:

for m from 0 to k0-1 do
gu:=gu+subs({n=n0,k=m},F):
od:
gu:

end:


#yakhas(hg,n): the asymptoic ratio of hg(n+1)/hg(n)
yakhas:=proc(hg,n) local lu,lu1,lu2,k:
lu:=normal(simplify(subs(n=n+1,hg)/hg)):
lu1:=expand(numer(lu)):
lu2:=expand(denom(lu)):
if degree(lu1,n)>degree(lu2,n) then
 RETURN(infinity):
fi:

if degree(lu1,n)<degree(lu2,n) then
 RETURN(0):
fi:
k:=degree(lu1):
coeff(lu1,n,k)/coeff(lu2,n,k):

end:


#Constant1(R1,R2,n,N): Given an acclerating pair
#finds the constant up to error 1/10^N in the 
#series term
Constant1:=proc(R1,R2,n,N)
local i,lu,T,mu:
if abs(yakhas(R1,n))>abs(yakhas(R2,n)) then
 T:=R2:
 mu:=0:
else
 T:=R1:
 mu:=evalf(subs(n=0,R1)):
fi:
lu:=1:
for i from 1 while abs(lu)>1/10^(N+5) do
lu:=evalf(subs(n=i,T)):
mu:=mu+lu:
od:

mu:

end:



#WZkummer(n,k,b,bitui): A WZ pair obtained from Kummer's 2F1(-1)
#identity, obtained by subsitituting to n bitui and  b is a number
#or paremeter, For example, try: WZkummer(n,k,0,n); 
WZkummer:=proc(n,k,b,bitui)
local F:
F:=subs(n=bitui,(-1)^k*k!*(k-b)!/(2*n+k+1)!/(2*n+k+1+b)!):

WZpair(F,k,n):
end:




#WZdixon(n,k,a,c,bitui): The generalizations of the WZ pair for Zeta(3)
#obtained by shadowing the WZ-pair implied by Dixon's identity
#with subsitituting to k bitui (starting with G)
#For example, try: WZdixon(n,k,0,0,k)
WZdixon:=proc(n,k,a,c,bitui)
local G,lu:
G:=subs(k=bitui,(n-k)!*(n+a)!*(n+c)!/(n+a+k+1)!/(n+1)!/(n+a-c+1)!):
lu:=WZpair(G,n,k):
lu[2],lu[1]:
end: