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

read EKHAD:

Digits:=100:

print(`Expanded package: Aug. 2021 `):

print(`This is AperyWZnew.txt, one of the packages that accompany the article`):
print(` Experimenting with Apery Limits and WZ pairs `):
print(``):
print(`By Robert Dogherty-Bliss  and Doron Zeilberger`):
print(``):

print(` Extending the of Original package : June 13, 2002`):
print(`First version, ca. Sept. 12, 2001`):

print(`This is AperyWZnew.txt, a Maple package that extends the work 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, with additions by Robert Dogherty-Bliss`):
print(``):
print(`For general help type: Ezra();`):
 print(``):

print(`For supporting procedures type: Ezra1();`):

print(`For the new procedures type: EzraN();`):

print(`For  procedures about specific WZ pairs type: EzraWZ();`):


EzraWZ:=proc()
print(`Theprocedures with WZ pairs, type: `):
print(` GWZzeta2, WZlog1x, WZlog2, WZzeta2, WZzeta3, WZchu, WZdixon, WZkummer, WZkummer0`):

end:



Ezra1:=proc()
print(`The supporting procedures are `):
print(`MulOpe, SeqFromRec `):

end:



EzraN:=proc()
print(`The New procedures are `):
print(`CSconj, zeilWZPnew, zeilWZPnewV  `):

end:

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,  MulOpe, Natcn`):
print(` WZpair, yakhas, zeilWZP, zeilWZPv `):

fi:

if nops([args])=1 and op(1,[args])=Acc then
print(`Acc(P,k,n): the acceleration formula implied by the WZ`):
print(`pair P=[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(P,bnk,n,k,n0): Given a WZ pair: P=[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(k,n),binomial(n,k)*binomial(n+k,k),n,k,5);`):
print(`try an(WZzeta2(k,n),binomial(n,k)^2*binomial(n+k,k),n,k,5);`):
print(`try an(WZzeta3(k,n),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(P,k,n): checks whether the pair P=[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(P,bnk,n,k,n0): Given a WZ-pair, P=[F,G], and  `):
print( ` averaging function, bnk (all expressions in n and k) `):
print(` and an integer  n0, finds the rational number `):
print(` an(P,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(k,n), binomial(n,k)^2*binomial(n+k,k),n,k,5);`):
fi:

if nops([args])=1 and op(1,[args])=cnk then
print(`cnk(P,n,k,n0,k0): Given a WZ pair P=[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(k,n),n,k,5,3);`):
fi:

if nops([args])=1 and op(1,[args])=cnkS then
print(`cnkS(P,n,k,m): Given a WZ pair P=[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(k,n),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])=CSconj then
print(`CSconj(d,n,N): Proof of a weaker version of the Chamberland-Straub conj. 9 in their nice article`):
print(` "Apery Limits Experiments and Proofs", arxXiv:2011.03400v1 6Nov 2020 for the`):
print(`sum of the d-th power of the n-th row of Pascal's triangle. d has to be between 3 and 9.`):
print(`it also states a technical lemma that would prove the whole conjecture (but we can't prove it yet.`) :
print(`it returns the recurrence operator in the variable n, and shift operator N`):
print(`and intial conditions whose proof of subdominance would imply the fully conjecture for that d.Try:`):
print(`CSconj(3,n,N):`):
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(P,b,n,k,K,N): Given a WZ pair P=[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(P,bnk,n,k,n0), of approximants to K`):
print(`according to Apery's method`):
print(` For example try: `):
print(`deltaWZb(WZlog2(k,n),binomial(n,k)*binomial(n+k,k),n,k, `):
print( `evalf(log(2)),30);` ):
print(`deltaWZb(WZzeta2(k,n),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, evalf(Zeta(3)),30); `):
fi:

if nops([args])=1 and op(1,[args])=GWZzeta2 then
print(`GWZzeta2(k,n,a,b,c): 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)), a,b,c are the parameters`):
print(`Try: `):
print(`GWZzeta2(k,n,0,0,0);`):
fi:


if nops([args])=1 and op(1,[args])=MulOpe then
print(`MulOpe(ope1,ope2,n,N): inputs two linear recurrence operators ope1(n,N), ope2(n,N), outputs the product ope1(n,N)ope2(n,N). Try:`):
print(`MulOpe(N+1,n/N+1,n,N);`):
fi:


if nops([args])=1 and op(1,[args])=Natcn then
print(` Natcn(P,n,k,n0): Given a WZ pair: P=[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(k,n),n,k,5); `):
fi:

if nops([args])=1 and op(1,[args])=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): inputs an operator given in temrs of`):
print(`N, and a list Ini, givining the values at n=1..,d`):
print(`to n=K. Try:`):
print(`SeqFromRec(N^2-N-1,n,N,[1,1],20);`):
fi:


if nops([args])=1 and op(1,[args])=WZchu then
print(`WZchu(k,n,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(k,n,0,0,0,n);`):
fi:


if nops([args])=1 and op(1,[args])=WZdixon then
print(`WZdixon(k,n,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(k,n,0,0,k);`):
fi:

if nops([args])=1 and op(1,[args])=WZkummer then
print(`WZkummer(k,n,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(k,n,0,n); `):
fi:


if nops([args])=1 and op(1,[args])=WZlog1x then
print(`WZlog1x(k,n,x)): The generalizations of the WZ pair`):
print(` for log((1+x)/x)`):
print(`For example, try: WZlog1x(k,n,x);`):
fi:


if nops([args])=1 and op(1,[args])=WZlog2 then
print(`WZlog2(k,n,x)): The generalizations of the WZ pair`):
print(` for log(2)`):
print(`For example, try: WZlog2(k,n);`):
fi:

if nops([args])=1 and op(1,[args])=WZzeta2 then
print(`WZzeta2(k,n)): The  WZ pair`):
print(` for Apery's proof for Zeta(2). Try:`):
print(` WZzeta2(k,n);`):
fi:

if nops([args])=1 and op(1,[args])=WZzeta3 then
print(`WZzeta3(k,n)): The  WZ pair`):
print(` for Apery's proof for Zeta(3). Try:`):
print(` WZzeta3(k,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,P,k,n,N), given an averaging function bnk`):
print(`and a WZ pair, P=[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),k,n,N); `):
print(` zeilWZP(binomial(n,k)^2*binomial(n+k,k),WZzeta2(k,n),k,n,N); `):
print(` zeilWZP(binomial(n,k)^2*binomial(n+k,k)^2,WZzeta3(k,n),k,n,N); `):

fi:

 
if nops([args])=1 and op(1,[args])=zeilWZPnew then
print(`zeilWZPnew(bnk,P,k,n,N,K), `):
print(`inputs:`):
print(` (i) an averaging function bnk `):
print(` (ii) P=[F,G], a WZ pair (all expressions in n,k) and a shift `):
print(` (iii) variable k `):
print(` (iv): variable n `):
print(` (v): Shift operator in n called N `):
print(` (vi): A large positive integer K (say 1000) `):
print(` Outputs a list whose entries are `):
print(` (i) The constant in question, Sum(F(n,0),n=0..infinity) `):
print(` (ii) The operator ope (in n and N) annihilating BOTH numerator and denominator `):
print(` (it may not be the minimal order operator for the denominator) `):
print(` (iii) the initial conditions for the numerator and denominator `):
print(` (iv) The smallest empirical delta of a(n)/b(n) for n=K-10..K `):
print(` Try: `):
print(` zeilWZPnew(binomial(n,k)*binomial(n+k,k),WZlog2(k,n),k,n,N,300); `):
print(` zeilWZPnew(binomial(n,k)^2*binomial(n+k,k),WZzeta2(k,n),k,n,N,300); `):
print(` zeilWZPnew(binomial(n,k)^2*binomial(n+k,k),WZchu(k,n,1/2,1/2,0,2*n),k,n,N,300); `):
print(` zeilWZPnew(binomial(n,k)^2*binomial(n+k,k)^2,WZzeta3(k,n),k,n,N,300); `):
print(` zeilWZPnew(binomial(n,k)^2*binomial(n+k,k),WZchu(k,n,1/2,1/2,0,n),k,n,N,300); `):
fi:

 
if nops([args])=1 and op(1,[args])=zeilWZPnewV then
print(`zeilWZPnewV(bnk,P,k,n,N,K), `):
print(`Verbose version of zeilWZPnew(bnk,P,k,n,N,K) (q.v.) `):
print(` Try: `):
print(` zeilWZPnewV(binomial(n,k)*binomial(n+k,k),WZlog2(k,n),k,n,N,300); `):
print(` zeilWZPnewV(binomial(n,k)^2*binomial(n+k,k),WZzeta2(k,n),k,n,N,300); `):
print(` zeilWZPnewV(binomial(n,k)^2*binomial(n+k,k)^2,WZzeta3(k,n),k,n,N,300); `):
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),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:=proc(k,n,x):
[(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)]:
end:

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

WZzeta2:=proc(k,n):
simplify([(-1)^(n+k)*k!^2*(n-k-1)!/(n+k+1)!,(-1)^(n+k)*k!^2*(n-k)!/(n+k+1)!*2/(n+1)]):
end:

WZzeta3:=proc(k,n):
simplify([(-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]):
end:

WZkummer0:=proc(k,n,b):
[
(-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)
]:
end:

###

#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(P,k,n) local F,G:
F:=P[1]: G:=P[2]:
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(P,n): the acceleration formula implied by the WZ
#pair P=[F,G] in the variables k,n
Acc:=proc(P,k,n) local F,G:
F:=P[1]: G:=P[2]:
if CheckWZ(P,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(P,n,k,m) local F,G:
F:=P[1]:
G:=P[2]:
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(P,bnk,n,k,n0): Given a WZ pair: P=[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(P,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(P,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(P,bnk,n,k,n0): Given a WZ-pair,P=[F,G] and 
#averaging function, bnk (all expressions in n and k)
#and an integer  n0, finds the rational number
#an(P,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(P,bnk,n,k,n0):an(P,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(P,b,n,k,K,N)
local lu,n0:
lu:=[seq(cn(P,b,n,k,n0),n0=1..N)]:
delta(lu,K):
end:


#Natcn(F,G,n,k,n0): Given a WZ pair: P=[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(P,n,k,n0)
local gu1,gu2,k0,mu:
gu1:=0: gu2:=0:
for k0 from 0 to n0 do
 mu:=cnk(P,n,k,n0,k0):
 gu1:=gu1+mu*denom(mu):
 gu2:=gu2+denom(mu):
od:
gu1/gu2:
end:

 
#zeilWZPv(bnk,P,k,n,N), given an averaging function bnk
#and a WZ pair,P= [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,P,k,n,N)
local F,G,i,j,gu,lu,R,cert,ku:

F:=P[1]:
G:=P[2]:
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,P,k,n,N), given an averaging function bnk
#and a WZ pair, P=[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,P,k,n,N)
local i,j,gu,lu,R,cert,ku,F,G:
F:=P[1]:
G:=P[2]:
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(P, `is not a WZ pair `):
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(k,n,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(k,n,a,b,c)
local F:
F:=(-1)^k*k!*(k+c)!*(n-k-1+a)!/(n+k+1+b)!:
simplify(WZpair(F,k,n)):
end:

#WZchu(k,n,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(k,n,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(P,n,k,n0,k0) local F,G,m,gu:
F:=P[1]:
G:=P[2]:

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:
eval(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(k,n,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(k,n,b,bitui)
local F:
F:=(-1)^k*k!*(k-b)!/(2*n+k+1)!/(2*n+k+1+b)!:

#F:=subs(k=bitui,(-1)^k*k!*(k-b)!/(2*n+k+1)!/(2*n+k+1+b)!):

WZpair(F,k,n):
end:




#WZdixon(k,n,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(k,n,0,0,k)
WZdixon:=proc(k,n,a,c,bitui)
local G,lu:
G:=(-n+k)!*(k+a)!*(k+c)!/(n+a+k+1)!/(k+1)!/(a-c+k+1)!:
G:=subs(k=bitui,G):
WZpair(G,k,n):

end:



#MulOpe(ope1,ope2,n,N): inputs two linear recurrence operators ope1(n,N), ope2(n,N), outputs the product ope1(n,N)ope2(n,N). Try:
#MulOpe(N+1,n/N+1,n,N);
MulOpe:=proc(ope1,ope2,n,N) local ope,i:
ope:=0:

for i from ldegree(ope1,N) to degree(ope1,N) do
 ope:=expand(ope+coeff(ope1,N,i)*subs(n=n+i,ope2)*N^i):
od:

add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):
end:



#SeqFromRec(ope,n,N,Ini,K): inputs an operator given in temrs of
#N, and a list Ini, givining the values at n=1..,d
#to n=K. Try:
#SeqFromRec(,n,N,[0,1],0,20);
SeqFromRec:=proc(ope,n,N,Ini,K) local d,i,T,n1,ope1:


d:=degree(ope,N):

ope1:=expand(subs(n=n-d,ope)/N^d):
ope1:=expand(ope1/coeff(ope1,N,0)):

ope1:=1-ope1:
ope1:=add(factor(coeff(ope1,N,-i))*N^(-i),i=1..d):

for n1 from 1 to nops(Ini) do
 T[n1]:=Ini[n1]:
od:

for n1 from 1+nops(Ini) to K do
 T[n1]:=add( subs(n=n1,coeff(ope1,N,-i))*T[n1-i],i=1..d):
od:

[seq(T[n1],n1=1..K)]:

end:


 
#zeilWZPnew(bnk,F,G,k,n,N,K), 
#inputs:
#(i) an averaging function bnk
#(ii) F,G, such that (F,G) is a WZ pair (all expressions in n,k) and a shift
#(iii) variable k
#(iv): variable n
#(v): Shift operator in n called N
#(vi): A large positive integer K (say 1000)
#Outputs a list whose entries are
#(i) The constant in question, Sum(F(n,0),n=0..infinity)
#(ii) The operator ope (in n and N) annihilating BOTH numerator and denominator
#(it may not be the minimal order operator for the denominator)
#(iii) the initial conditions for the numerator and denominator
#(iv) The smallest empirical delta of a(n)/b(n) for n=K-10..K
#Try:
#zeilWZPnew(binomial(n,k)*binomial(n+k,k),WZlog2(k,n),k,n,N,100);

zeilWZPnew:=proc(bnk,P,k,n,N,K) local mu,gu,ope,Ini1,Ini2,C,SeqA,SeqB,n0,d,del,S,i:
Digits:=1000:

 gu:=zeilWZP(bnk,P,k,n,N):

mu:=Acc(P,k,n):

C:=sum(mu[1],n=0..infinity):

ope:=MulOpe(gu[3],gu[1],n,N):

d:=degree(ope,N):



Ini1:=simplify([seq(an(P,bnk,n,k,n0),n0=1..d)]):

Ini2:=[seq(bn(bnk,n,k,n0),n0=1..d)]:

if not (type(Ini1[1],integer) or type(Ini1[1],fraction)) then
RETURN([C,ope,Ini1,Ini2]):

fi:



SeqA:=SeqFromRec(ope,n,N,Ini1,K):
SeqB:=SeqFromRec(ope,n,N,Ini2,K):



S:=[seq(SeqA[i]/SeqB[i],i=K-10..K)]:


del:=min(seq(evalf(-log(abs(evalf(C)-S[i]))/log(denom(S[i])))-1,i=1..nops(S))):
del:=evalf(del,10):

[C,ope,Ini1,Ini2,del]:

end:



#zeilWZPnewV(bnk,P,k,n,N,K): verbose version of zeilWZPnew(bnk,P,k,n,N,K) (q.v.)
zeilWZPnewV:=proc(bnk,P,k,n,N,K) local gu,ope,Ini1,Ini2,del,x,a,b,i:

gu:=zeilWZPnew(bnk,P,k,n,N,K):

ope:=gu[2]:

Ini1:=gu[3]:

Ini2:=gu[4]:

del:=gu[5]:

if nops(gu)<5 then
 RETURN(FAIL):
fi:
print(``):

print(`Diophantine Approximations to the Constant`, gu[1], `in terms of ratios of sequences satisfying the same recurrence `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let a(n), b(n) be the sequences satisfying the linear recurrence `):

print(``):
print(add(coeff(ope,N,i)*x(n+i),i=0..degree(ope,N))=0):
print(``):

print(`and in Maple format`):

print(``):
lprint(add(coeff(ope,N,i)*x(n+i),i=0..degree(ope,N))=0):
print(``):

print(`subject to the initial conditions `):

print(``):
print(seq(a(i)=Ini1[i],i=1..nops(Ini1))):
print(``):

print(``):
print(seq(b(i)=Ini2[i],i=1..nops(Ini2))):
print(``):
print(`Then the limit a(n)/b(n), as n goes to infinity is indeed the constant`, gu[1]):
print(``):
print(`The delta such that the error is less than the reciporocal of the denominator to the power delta+1 is approximately`, del):
print(``):

if del>0 then
 print(`This indicated a proof that`, gu[1], `is irrational with an irrationality measure about`, evalf(1+1/del,10) ):
elif
del>-0.7 then
print(`This is nice, but not good enough to prove irrationality`):
else
print(`This is very disappointing`):
fi:

end:


#CSconj(d): Proof of a weaker version of the Chamberland-Straub conj. 9 in their nice article
# "Apery Limits Experiments and Proofs", arxXiv:2011.03400v1 6Nov 2020 for the
#sum of the d-th power of the n-th row of Pascal's triangle. try:
#CSconj(3);

CSconj:=proc(d,n,N) local ope,ord,INI1,INI1a, T,n1,ope1,i,k,gu,lu,INI2,INI2a,INI3,opeK,delta,E,INI4,CONST,m,hef,Bseq:

if not (type(d,integer) and d>=3) then

 print( d, `must be an integer at least 3`):

 RETURN(FAIL):
fi:

if d>=10 then
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:

ope:=zeil(binomial(n,k)^d,k,n,N)[1]:

ord:=degree(ope,N):

ope1:=ope/coeff(ope,N,ord):

ope1:=1-expand(subs(n=n-ord,ope1)/N^ord):

ope1:=add(factor(coeff(ope1,N,-i))*N^(-i),i=1..ord):

for n1 from -(ord-2) to 0 do
T[n1]:=0:
od:

T[1]:=1:

for n1 from 2 to ord do
 T[n1]:=add(subs(n=n1,coeff(ope1,N,-i))*T[n1-i],i=1..ord):
od:

INI1:=[seq(T[n1],n1=1..ord)]:
INI2:=[seq(add(binomial(n1,k)^d,k=0..ord),n1=1..ord)]:



gu:=zeilWZPnew(binomial(n,k)^d,WZzeta2(k,n)/(d+1),k,n,N,1000):

opeK:=gu[2]:
INI3:=gu[3]:


print(`A recurrence satisfied by the Sum of the`, d, `-th power of the n-th row of the Pascal triangle and another solution`):
print(`of that same recurrence, such that the ratios converge to`, Zeta(2)/(d+1) ):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Fact: Let `):

print(B(n)=Sum(binomial(n,k)^d,k=0..n)):
print(``):
lu:=cnkS(WZzeta2(k,n)/(d+1),n,k,m):
print(A(n)=Sum(binomial(n,k)^d*lu ,k=0..n)):
print(``):
print(`then A(n)/B(n) converges to`, Zeta(2)/(d+1) ):
print(``):
print(`The potential function `):
print(``):
print(lu):
print(``):
print(`obviously converges to`, Zeta(2)/(d+1), ` as (n,k) go to infinity (in each direction) and A(n)/B(n) are weighted averages hence must also`):
print(``):
print(`converge to that limit `):
print(``):
print(`Using the Zeilberger algorithm and its extension it follows that both B(n) and A(n) satisfy the SAME linear recurrence `):
print(``):
print(add(coeff(opeK,N,i)*X(n+i),i=0..degree(opeK,N))=0):
print(``):

INI2a:=SeqFromRec(ope,n,N,INI2,degree(opeK,N)):
INI1a:=SeqFromRec(ope,n,N,INI1,degree(opeK,N)):



print(``):
print(`with the initial conditions, respectively `):
print(``):
print(seq(B(i)=INI2a[i],i=1..degree(opeK,N))):
print(``):
print(seq(A(i)=INI3[i],i=1..degree(opeK,N))):
print(``):
print(`Hence we found a linear recurrence (albeit not the minimal one) satisfied by the B(n), and another solution A(n) of that`):
print(`very same recurrence, with different initial conditions such that their ratio converges to`, Zeta(2)/(d+1) ):
print(``):
print(`Furthemore `):
print(``):
print(abs(A(n)/B(n)- Zeta(2)/(d+1))<=CONST/n^(delta)):
print(``):
print(`where delta is about`, 1+gu[5] ):
print(``):
print(`This proves a weaker form of a conjecture of Marc Chamberland and  Armin Straub`):
print(` "Apery Limits Experiments and Proofs", arxXiv:2011.03400v1 6 Nov 2020.`):
print(``):
print(`They conjectured that the sequence with initial conditions`):
print(``):
print(seq(C(n1)=0,n1=-(ord-2)..0), C(1)=1 ):
print(``):
print(`satisfying the MINIMAL recurrence satisfied by B(n), namely `):
print(``):
print(add(coeff(ope,N,i)*X(n+i),i=0..degree(ope,N))=0):
print(``):
print(`also has the property that the limit of C(n)/B(n) converges to`, Zeta(2)/(d+1)):
print(``):
print(`It is readily seen that C(n) also satisfies the  (non-minimal) recurrence above, and hence the difference`):
print(E(n)=C(n)-A(n)):
print(``):
INI4:=INI1a-INI3:
print(`satisfied that recurrence, with initial conditions`):
print(``):
print(seq(E(n1)=INI4[n1],n1=1..degree(opeK,N)) ):
print(``):
print(`The full conjecture of Chamberland and Straub would follow once that we can prove that this sequence of differences is SUB-DOMINANT`):
print(``):
print(`in other words that the limit of `, E(n)/B(n), `goes to 0 as n goes to infinity `):
print(``):
print(`We confess that we do not know how to prove this rigorously, but let's do it empirically `):
print(``):
hef:=SeqFromRec(opeK,n,N,INI4,1000):
Bseq:=SeqFromRec(ope,n,N,INI2,1000):

print(``):
print(`The terms from n=991 to n=1000 are `):
print(``):
print(seq(evalf(hef[i]/Bseq[i]),i=990..1000)):
print(``):

[opeK,INI4]:

end: