#####################################################################
## Zudilin.txt Save this file as  Zudilin.txt to use it,             #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  `Zudilin.txt`: <Enter>                      #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Robert Dougherty-Bliss and                             # 
#Doron Zeilberger, Rutgers University ,                              #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(` Written:  May 2022: tested for Maple 2020 `):
print(`This is Zudilin.txt, A Maple package`):
print(`accompanying Robert Dougherty-Bliss and Doron Zeilberger's  article: `):
print(` Exploring General Apery Limits via the Zudilin-Straub t-transform  `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Zudilin.txt  .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print(`----------------------------------------`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(``):
 print(`----------------------------------------`):
print(``):

print(`For a list of the supporting functions;`):
print(` type "ezra1();". For specific help type "ezra(procedure_name);" `):
print():

 print(`----------------------------------------`):
print(``):
print(`For a list of the General functions;`):
print(` type "ezraG();". For specific help type "ezra(procedure_name);" `):
print(``):

 print(`----------------------------------------`):


print(`For a list of the STORY functions;`):
print(` type "ezraSt();". For specific help type "ezra(procedure_name);" `):


 print(`----------------------------------------`):

Digits:=1000:


ezraSt:=proc()
if args=NULL then

print(`The STORY procedures are`):
print(`  Paper, Paper2, PaperG `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  ApplyOpe, cnk, cnkD, cnk2, cnkF, CrL, delt, EvalAjnk, Eval1RF, EvalRF, EvalZTF, FindMaxk,  rf1, SeqFromRec, SeqFromRecG, ZT1 `):

else
ezra(args):

fi:


end:



ezraG:=proc()
if args=NULL then

print(`The General procedures are`):
print(`   AZcG,   EvalFrb, zeilL1rG`):

else
ezra(args):

fi:


end:


ezra:=proc()
if args=NULL then

print(` Zudilin.txt: A Maple package for  exploring Apery Limits using the Zudilin-Straub  t-Transform`):

print(`The MAIN procedures are`):
print(` AZc, RFseq, RowZT, ZT, zeil, zeilL, zeilL1, zeilL1r, `):


elif nargs=1 and args[1]=ApplyOpe then
print(`ApplyOpe(ope,n,N,L): applies the operator ope in n, and N to the list L, Try:`):
print(`ApplyOpe(2-N,n,N,[seq(2^i,i=1..20)]);`):

elif nargs=1 and args[1]=AZc then
print(`AZc(F,k,n,r,M):  Inputs a `):
print(`hypergeometric term F (phrased in terms of the variables k and n, a small pos. integer r and a large pos. integer M.`):
print(`outputs the  approx. of the Apery-Zudilin `):
print(`constant of the coeff. of t^r the Zudilin-Straub transform of F, followed by the estimated delta.  Finally it reurns the number of digits in the approximaion with M terms. Try:`):
print(`AZc(binomial(n,k)^3,k,n,2,1000);`):



elif nargs=1 and args[1]=AZcG then
print(`AZcG(F,k,n,r,MaxC,M):  Inputs a `):
print(`hypergeometric term F (phrased in terms of the varlaibes k and n). It also inputs a small positive integer r, a medium-sized integer MaxC and`):
print(`a large integer M.`):
print(`It outputs the  approx. of the Apery-Zudilin `):
print(`constant of the coeff. of t^r the Zudilin-Straub transform of F, followed by the estimated delta. Try:`):
print(`AZcG(binomial(n,k)*binomial(n+k,k),k,n,2,15,1000);`):


elif nargs=1 and args[1]=CrL then
print(`CrL(A,r,L): The symbolic expression for  cnk(n,k,r,L) (q.v.) in terms of A[j]'s where it stands for A[j](n,k). Try:`):
print(`CrL(A,4,L);`):

elif nargs=1 and args[1]=cnk then
print(`cnk(n,k,r,L): The coefficient of the Taylor expansion of t^r in (1/((1+t)*...*(1+t/k))*((1-t)*...*(1-t/(n-k)))^L. Try:`):
print(`[seq(cnk(10,k,2,3),k=0..10)];`):


elif nargs=1 and args[1]=cnkD then
print(`cnkD(n,k,r,L): The coefficient of the Taylor expansion of t^r in (1/((1+t)*...*(1+t/k))*((1-t)*...*(1-t/(n-k)))^L,DONE DIRECTLY, via CrL(A,r,L) (q.v.). Try:`):
print(`[seq(cnkD(10,k,2,3),k=0..10)];`):

elif nargs=1 and args[1]=cnkF then
print(`cnkF(F,n,k,n1,k1,r): The coefficient of the Taylor expansion of t^r in  `):
print(`ZT(F,k,t,RF)/F for numeric n=n1 and k=k1. Try:`):
print(`cnkF(binomial(n,k)^3,n,k,5,3,2);`):

elif nargs=1 and args[1]=cnk2 then
print(`cnk2(n,k,L): The coefficient of the Taylor expansion of t^2 in (1/((1+t)*...*(1+t/k))*((1-t)*...*(1-t/(n-k)))^L. Using the closed-form formula`):
print(`[seq(cnk2(10,k,3),k=0..10)];`):

elif nargs=1 and args[1]=delt then
print(`delt(a,c): Given a rational number a and a constant`):
print(`c, finds the delta such that |a-c|=1/denom(a)^(1+delta)`):
print(`For example, try delta(22/7,evalf(Pi)):`):

elif nargs=1 and args[1]=EvalAjnk then
print(`EvalAjnk(j,n,k): Evaluates A_j(n,k):=1/j*((-1)^(j-1)*Sum(1/i^j,i=1..k) - Sum(1/i^j,i=1..n-k) ). Try:`):
print(`EvalAjnk(5,20,5);`):

elif nargs=1 and args[1]=EvalFrb then
print(`EvalFrb(r,b,N): the numeric value of the atomic quantity F_{r,b}(N):=(-1)^r/r*Sum(b^r/j^r,j=1..N). Try:`):
print(`EvalFrb(2,3,10);`):

elif nargs=1 and args[1]=Eval1RF then
print(`Eval1RF(F,RF,n,k,n1,k1): Given an expression in F and k, featuring RF outputs its evaluation at n=n1,k=k1`):
print(`Try: `):
print(`Eval1RF(ZT(n!/k!/(n-k)!,k,t,RF),RF,n,k,5,3);`):

elif nargs=1 and args[1]=EvalRF then
print(`EvalRF(F,RF,n,k,n1): Given an expression in F and k, featuring RF, outputs Sum(F(k1),k1=0..n1):`):
print(`Try: `):
print(`EvalRF(ZT(n!/k!/(n-k)!,k,t,RF),RF,n,k,5);`):

elif nargs=1 and args[1]=EvalZTF then
print(`EvalZTF(F,k,n,k1,n1,r): Evaluates  the coeff. of t^r in the Taylor expansion of  ZT(F,k,n,RF)/F at k=1,n=n1. Try`):
print(`EvalZTF(binomial(n,k)^3,k,n,10,20,2);`):

elif nargs=1 and args[1]=FindMaxk then
print(`FindMaxk(F,n,k): Inputs a hypergeometric term in n and k, outputs the number, let's call it alpha,  such that k=alpha*n t makes F(n,k) as large as possible. Try:`):
print(`FindMaxk(binomial(n,k),n,k);`):


elif nargs=1 and args[1]=Paper then
print(`Paper(K1,K2,M): Inputs two positive integers K1 (at least 3) and K2, and a large integer M (for the input to AZc (q.v.),  and outputs a paper with all the Apery limits gotten from`):
print(`considering the coefficient of t^r, for r=1...L-1, in the Zudilin-Straub transform of binomial(n,k)^L*a^k for L from 3 to K1 and a from 1 to K2, and r from 2 to L-1.`):
print(`Try:`):
print(`Paper(4,2,2000);`):

elif nargs=1 and args[1]=Paper2 then
print(`Paper2(K1,K2,M): Inputs two positive integers K1 (at least 3) and K2, and a large integer M (for the input to AZc (q.v.),  and outputs a paper with all the Apery limits gotten from`):
print(`considering the coefficient of t^2 in the Zudilin-Straub transform of binomial(n,k)^L*a^k for L from 3 to K1 and a from 1 to K2.`):
print(`Try:`):
print(`Paper2(4,2,2000);`):

elif nargs=1 and args[1]=PaperG then
print(`PaperG(F,k,n,r,MaxC,M):Inputs a hypergeometric term F in discrete variables n and k, a small positive integer r, a medium integer MexC and a large M`):
print(`gives information about the Apery constant inspired by the coeff. of t^t in the Straub-Zudilin transform  if Sum(F,k=0..n). Try:`):
print(`PaperG(binomial(n,k)^3,k,n,4,10,2000);`):

elif nargs=1 and args[1]=rf1 then
print(`rf1(a,k): a(a+1)...(a+k-1)`):

elif nargs=1 and args[1]=RFseq then
print(`RFseq(F,RF,n,k,m): The list of the first m terms of EvalRF(F,RF,n,k,n1) for n1=1 to n1=m. Try:`):
print(`RFseq(ZT(n!/k!/(n-k)!,k,t,RF),RF,n,k,10);`):

elif nargs=1 and args[1]=RowZT then
print(`RowZT(F,k,n,n1,r): Inputs a hypergeometric term F phrased in terms of variables k and n, an integer n1, and a small integer r`):
print(`outputs the list of length n1+1 whose (k1+1)-th entry is the coeff. of t^r in  ZT(F,k,t,RF). Try:`):
print(`RowZT(binomial(n,k)^3,k,n,5,2); `):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):



elif nargs=1 and args[1]=SeqFromRecG then
print(`SeqFromRecG(ope,n,N,Ini,RHS1, K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)] and a holnomic description of the right side`):
print(`satisfied by the recurrence ope(n,N)f(n)=RHS1(n)`):
print(`extends it to the first K values. Try:`):
print(`SeqFromRecG(N-2,n,N,[1],[N-1,[4]],30);`):

elif nargs=1 and args[1]=ZT then
print(`ZT(F,k,t,RF): Given a product/quotients of terms of the form (an+bk+c)! and (an-bk+c)! converts each to  RF(1+bt,an+bk+c) . Try: `):
print(`ZT(n!/(k!*(n-k)!),k,t,RF);`):

elif nargs=1 and args[1]=ZT1 then
print(`ZT1(F,k,n,t,RF) : converts (an+b*k+c)! to RF(1+b*t,F). Try:`):
print(`ZT1((5*n-3*k+5)!,k,t,RF);`):

elif nargs=1 and args[1]=zeil then
print(`zeil(F,k,n,N): outputs a hypergemeoetric F in k,n, and outputs an operator ope, and a certificate R(n,k), such that if G(n,k):=R(n,k)*F(n,k), then`):
print(` ope(n,N) F(n,k)=G(n,k+1)-G(n,k). Try`):
print(`zeil(n!^3/k!^3/(n-k)!^3,k,n,N);`):


elif nargs=1 and args[1]=zeilL then
print(`zeilL(F,k,n,t,N,RF): inputs a proper hypergeometric term in k,n, i.e. a product/quotient of terms of the form (an+bk+c)! and a^k outputs`):
print(`the pair`):
print(`[ope(n,N), RHS(n,t)] such if F1:=ZT(F,k,n,t,RF) is the Zudilin-Straub transform of F (in terms of t and the RAISING FACTORIAL RF)`):
print(`that  ope(n,N) Sum(F1(t,k),k=0..n)=RHS(n,t). Try:`):
print(`zeilL(binomial(n,k)^3,k,n,t,N,RF);`):

elif nargs=1 and args[1]=zeilL1 then
print(`zeilL1(F,k,n,t,N,M): inputs a proper hypergeometric term in k,n, i.e. a product/quotient of terms of the form (an+bk+c)! and a^k outputs,`):
print(` and a positive integer M, outputs the triple`):
print(`[ope(n,N), L0,L1, L2] where ope(n,N) is the recurrence operator gotten by the Zeilberger algorithm, L0 is the original (t=0) case, L1 are the first M terms of the sum of the Zudilin-Straub transform`):
print(`sum (using t), and L2 are the first M terms of the right side. Try:`):
print(`zeilL1(binomial(n,k)^3,k,n,t,N,10);`):

elif nargs=1 and args[1]=zeilL1r then
print(`zeilL1r(F,k,n,N,M,r): Same input as zeilL1(F,k,n,t,N,M) but instead of L2, L3, outputs the`):
print(`NUMERICAL sequence of the Taylor coefficients of t^r on both sides.  It also returns the quotient with the underlying binomial sum.Try:`):
print(`zeilL1r(binomial(n,k)^3,k,n,N,10,2);`):

elif nargs=1 and args[1]=zeilL1rG then
print(`zeilL1rG(F,k,n,N,MaxC,r): Like zeilL1r(F,k,n,N,M,r), but also tries to guess a holonomic deserciption to the right side. Try:`):
print(`zeilL1rG(binomial(n,k)^3*2^k,k,n,N,10,3);`):

else
 print(`There is no such thing as`, args):

fi:


end:




#Start FROM FindRec.txt

#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:



if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:

#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
# print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:


#End FROM FindRec.txt


#delt(a,c): Given a rational number a and a constant
#c, finds the delta such that |a-c|=1/denom(a)^(1+delta)
#For example, try delta(22/7,evalf(Pi)):
delt:=proc(a,c): 

if a=0 then
RETURN(FAIL):
fi:

evalf(-log(abs(c-a))/log(denom(a)))-1: end:

#Start from FindRec



#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:

ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:

#end from FindRec

####Start From EKHAD


with(SolveTools):
solve1:=Linear: 


rf:=proc(a,b):
 (a+b-1)!/(a-1)!:
end:


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

findgQR:=proc(Q,R,k,L)
local j,g:
 
for j from 1 to L do
 g:=gcd(Q,subs(k=k+j,R)):
 if degree(g,k)>0 then
  RETURN(g,j):
 fi:
od:
1,0:
end:



  
simplify1:=proc(bitu,k,a)
local gu,gu1,i,khez,sp:
sp:=1:
gu:=bitu:
 
if type(gu,`*`) then
 
 for i from 1 to nops(gu) do
  gu1:=op(i,gu):
 
 if type(gu1,`^`) and type(op(2,gu1), integer) then
 
     khez:=op(2,gu1):
     gu1:=op(1,gu1):
     sp:=sp*simplify(subs(k=k+a,gu1)/gu1)^khez:
  
   else 
     sp:=sp*simplify(subs(k=k+a,gu1)/gu1):
 
   fi:
 
od:
 
elif type(gu,`^`) and type(op(2,gu), integer) then
   khez:=op(2,gu):
   gu1:=op(1,gu):
   sp:=sp*simplify(subs(k=k+a,gu1)/gu1)^khez:  
 
 
else
 
 sp:=simplify(subs(k=k+a,gu)/gu):
 
fi:
 
sp:
 
end:

 
hafel:=proc(POL,SUMMAND,ope,n,N)
local i,FAC,degg,bi,rat,POL1,SUMMAND1,zub:
degg:=degree(ope,N):
FAC:=0:
 
for i from 0 to degg do
 bi:=coeff(ope,N,i):
 rat:=simplify1(SUMMAND,n,i):
 FAC:=FAC+bi*subs(n=n+i,POL)*rat:
 FAC:=normal(FAC):
od:
 
POL1:=numer(FAC):
zub:=denom(FAC):
SUMMAND1:=SUMMAND/zub:
RETURN(POL1,SUMMAND1,zub):
end:
 
 
ct:=proc(SUMMAND,ORDER,k,n,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,res1,hakhi,g,
A1, B1,  P1, SDZ, SDZ1,
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,ope1,denFAC,ope1a,apc,bpc:

 
if nargs<>5 then
ERROR(`Syntax: ct(SUMMAND,ORDER,summation_variable,auxiliary_var, Shift_n)`):
fi:
 
if tovu(convert(SUMMAND,factorial),k,n)=0 then
ERROR(`The summand can be separated into f(`,k,`)g(`,n,`)`): 
fi:
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od:
 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
denFAC:=gu[3]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
 
 
res1:=findgQR(Q,R,k,100):
 
while res1[2]<>0 do
 g:=res1[1]:
 hakhi:=res1[2]:
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=findgQR(Q,R,k,100):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
 
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,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:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
fi:
gu:={}:
for i1 from 0 to k1 do
gu:=gu union {apc[i1]=1}:
od:
for i1 from 0 to ORDER do
gu:=gu union {bpc[i1]=1}:
od:
fu:=subs(gu,fu):
ope:=subs(gu,ope):
 
ope:=pashet(ope,N):
ope1:=ope*denom(ope):
 
SDZ:=denom(ope)*subs(k=k+1,Q)*fu/P1/denFAC :
SDZ1:=subs(k=k-1,SDZ)*simplify1(convert(SUMMAND,factorial),k,-1):
SDZ1:=simplify(SDZ1):
 
ope1a:=0:
for i from 0 to degree(ope1,N) do
ope1a:=ope1a+sort(coeff(ope1,N,i)*N^i):
od:
 
RETURN(ope1a,SDZ1):
end:

 
 
zeil4:=proc(SUMMAND,k,n,N)
local ORDER,MAXORDER,gu,gu1,resh:
 
 
 MAXORDER:=20:
 
  resh:=[]:
 
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
 gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
 
 
zeil5:=proc(SUMMAND,k,n,N,MAXORDER)
local ORDER,gu,gu1,resh:
 
resh:=[]:
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
 gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
 
 
zeil6:=proc(SUMMAND,k,n,N,MAXORDER,resh)
local ORDER,gu,gu1:
 
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
 gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
zeil:=proc()
option remember: 
if nargs=4 then
 zeil4(args):
elif nargs=5 then
 zeil5(args):
elif nargs=6 then
 zeil6(args):
else
  print(`zeil(SUMMAND,k,n,N) or zeil(SUMMAND,k,n,N,MAXORDER) or`):
  ERROR(`zeil(SUMMAND,k,n,N,MAXORDER,param_list)`):
fi:
 
 
 
end:




 
tovu:=proc(SU,k,n)
local gu:
 
gu:=simplify1(simplify1(SU,n,1),k,1):
 
if gu=1 then
RETURN(0):
else
RETURN(1):
fi:
end:


 
cttry:=proc(SUMMAND,ORDER,k,n,resh,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,res1,hakhi,g,
A1, B1,  P1, 
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,eqg, apc, bpc:
 
if nargs<>6 then
ERROR(`The syntax is cttry(term,ORDER,sum_variable,aux_var,para_list,Shift)`):
fi:
 
if tovu(convert(SUMMAND,factorial),k,n)=0 then
ERROR(`The summand can be separated into f(`,k,`)g(`,n,`)`): 
fi:
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od:
 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
 
res1:=findgQR(Q,R,k,100):
 
while res1[2]<>0 do
 g:=res1[1]:
 hakhi:=res1[2]:
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=findgQR(Q,R,k,100):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,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:
eqg:=subs(n=1/2,eq):
 
for i from 1 to nops(resh) do
 eqg:=subs(op(i,resh)=i^2+3,eqg):
od:
 
va1:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
 else
  RETURN(1):
fi:
 
end:


###End From EKHAD 

#ZT1: converts (an+b*k+c)! to RF(1+b*t,F). Try:
#ZT1((5*n-3*k+5)!,k,t,RF);

ZT1:=proc(F,k,t,RF)

if type(F,numeric) or (type(F,`^`) and not type(op(1,F),function) ) then
 RETURN(F):
fi:

if type(F,function) and op(0,F)=factorial then
RETURN(RF(1+coeff(op(1,F),k)*t,op(1,F))):
fi:


if type(F,`^`) and type(op(1,F),function) and op(0,op(1,F))=factorial then
RETURN( RF(1+coeff(op(1,op(1,F)),k)*t,op(1,op(1,F)))^op(2,F)):
fi:

FAIL:

end:

#ZT(F,k,t,RF): Given a product/quotients of terms of the form (an+bk+c)! and (an-bk+c)! converts each to  RF(1+bt,an+bk+c) . Try
#ZT(n!/(k!*(n-k)!),k,t,RF);
ZT:=proc(F,k,t,RF) local i,F1:

F1:=convert(F,factorial):

if not type(F1,`*`) then
 RETURN(ZT1(F1,k,t,RF)):
else
 RETURN(mul(ZT1(op(i,F1),k,t,RF),i=1..nops(F1))):
fi:

end:


#rf1(a,k): a(a+1)...(a+k-1)
rf1:=proc(a,k) local i: mul(a+i,i=0..k-1):end:


#Eval1RF(F,RF,n,k,n1,k1): Given an expression in F and k, featuring RF, outputs its evaluation with n=n1,k=k1
#Try
#Eval1RF(ZT(n!/k!/(n-k)!,k,t,RF),RF,n,k,5,3);
Eval1RF:=proc(F,RF,n,k,n1,k1):

normal(eval(subs(RF=rf1,subs({n=n1,k=k1},F)))):

end:


EvalRF:=proc(F,RF,n,k,n1) local k1:
normal(add(Eval1RF(F,RF,n,k,n1,k1),k1=0..n1)):
end:


#RFseq(F,RF,n,k,m): The list of the first m terms of EvalRF(F,RF,n,k,n1) for n1=1 to n1=m. Try
#RFseq(ZT(n!/k!/(n-k)!,k,t,RF),RF,n,k,10);
RFseq:=proc(F,RF,n,k,m) local n1:
normal([seq(EvalRF(F,RF,n,k,n1),n1=1..m)]):
end:



 

#ApplyOpe(ope,n,N,L): applies the operator ope in n, and N to the list L, Try:
#ApplyOpe(2-N,n,N,[seq(2^i,i=1..20)]);
ApplyOpe:=proc(ope,n,N,L) local i,n1:
normal([seq(add(subs(n=n1,coeff(ope,N,i))*L[n1+i],i=0..degree(ope,N)),n1=1..nops(L)-degree(ope,N))]):
end:



#zeilL(F,k,n,t,N,RF): inputs a proper hypergeometric term in k,n, i.e. a product/quotient of terms of the form (an+bk+c)! and a^k outputs
#the pair
#[ope(n,N), RHS(n,t)] such if F1:=ZT(F,k,n,t,RF) is the Zudilin-Straub transform of F (in terms of t and the RAISING FACTORUAL RF)
#that  ope(n,N) Sum(F1(t,k),k=0..n)=RHS(n,t). Try:
#zeilL(binomial(n,k)^3,k,n,t,RF);
zeilL:=proc(F,k,n,t,N,RF) local F1,ope,YEMIN,G:
F1:=ZT(F,k,t,RF):
ope:=zeil(eval(subs(RF=rf,F1)),k,n,N):
G:=ope[2]:
ope:=ope[1]:

YEMIN:=subs(k=n+1,G)*subs(k=n+1,F1)-subs(k=0,G)*subs(k=0,F1):

ope, YEMIN:

end:




#zeilL1(F,k,n,t,N,M): inputs a proper hypergeometric term in k,n, i.e. a product/quotient of terms of the form (an+bk+c)! and a^k outputs#
#and a positive integer M, outputs the triple
#[ope(n,N), L1, L2] where ope(n,N) is the recurrence operator gotten by the Zeilberger algorithm, L1 are the first M terms of the sum of the Zudilin-Straub transform
#sum (using t), and L2 are the first M terms of the right side. Try:
#zeilL1(binomial(n,k)^3,k,n,N,RF,10);
zeilL1:=proc(F,k,n,t,N,M) local F1,ope,L1,L2,RF:
F1:=ZT(F,k,t,RF):
ope:=zeil(F,k,n,N)[1]:
L1:=RFseq(F1,RF,n,k,M):

if subs(t=0,L1)<>RFseq(F,RF,n,k,M) then
 RETURN(FAIL):
fi:

L2:=factor(ApplyOpe(ope,n,N,L1)):
L1:=factor(L1):

if subs(t=0,L2)<>[0$nops(L2)] then
 RETURN(FAIL):
fi:

[ope,subs(t=0,L1), L1,L2]:
end:





#zeilL1r(F,k,n,N,M,r): Same input as zeilL1(F,k,n,N,N) but instead of L2, L3, outputs the
#NUMERICAL sequence of the Taylor coefficients of t^r on both sides. It also returns the sequence of quotient with the unerlying binomial sum. Try:
#zeilL1r(binomial(n,k)^3,k,n,N,10,2);
zeilL1r:=proc(F,k,n,N,M,r) local t,gu,i,L1,L2,ope,L0:
option remember:
gu:=zeilL1(F,k,n,t,N,M):

if gu=FAIL then
  RETURN(FAIL):
fi:

ope:=gu[1]:


L0:=gu[2]:

L1:=gu[3]:

L2:=gu[4]:

L0:=subs(t=0,L1):

L1:= [seq(coeff(taylor(L1[i],t=0,r+6),t,r),i=1..nops(L1))]:
L2:= [seq(coeff(taylor(L2[i],t=0,r+6),t,r),i=1..nops(L2))]:

[ope,L0,L1,L2,[seq(evalf(L1[i]/L0[i],10),i=1..M)]]:

end:









#RowZT(F,k,n,n1,r): Inputs a hypergeometric term F phrased in terms of variables k and n, an integer n1, and a small integer r
#outputs the list of length n1+1 whose (k1+1)-th entry is the coeff. of t^r in 
#RowZT:=proc(F,k,n,n1,r)
RowZT:=proc(F,k,n,n1,r) local RF,t,F1,F1a,k1,gu:
F1:=ZT(F,k,t,RF):
gu:=[]:

for k1 from 0 to n1 do
 F1a:=subs({n=n1,k=k1},F1):
 F1a:=normal(eval(subs(RF=rf1,F1a))):
 F1a:=F1a/subs(t=0,F1a):
  F1a:=coeff(taylor(F1a,t=0,r+3),t,r):
 gu:=[op(gu),F1a]:
od:

gu:

end:




#AZc(F,k,n,r,M):  Inputs a 
#hypergemetric term F (phrased in terms of the varlaibes k and n. It also outputs the estimated delta. a small positive integer r and a large positive integr M,
#outputs the  approx. of the Apery-Zudilin 
#constant of the coeff. of t^r the Zudilin-Straub transform of F, followed by the estimated delta. Finally it reurns the number of digits in the approximaion with M terms. Try:
#AZc(binomial(n,k)^3,k,n,2,1000);
AZc:=proc(F,k,n,r,M) local gu,N,ope,L0,L1,C1,C2,D1,Ini0,Ini1,i,lu:

if M<100 then
 print(M, `should have been at least 100`):
 RETURN(FAIL):
fi:

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

gu:=zeilL1r(F,k,n,N,degree(ope,N)+10,r):

if gu=FAIL then
 RETURN(FAIL):
fi:

 if gu[4]<>[0$nops(gu[4])] then
# print(`Use AZcG`):
 RETURN(FAIL):
fi:

ope:=gu[1]:

L0:=gu[2]:

L1:=gu[3]:

Ini0:=[op(1..degree(ope,N),L0)]:
Ini1:=[op(1..degree(ope,N),gu[3])]:





L0:=SeqFromRec(ope,n,N,Ini0,M):
L1:=SeqFromRec(ope,n,N,Ini1,M):


C1:=L1[M]/L0[M]:
C2:=L1[M-1]/L0[M-1]:
 
if abs(C1-C2)>0 then

D1:=trunc(-log[10](abs(C1-C2))):
else
D1:=FAIL:

fi:


if D1<>FAIL and D1<20 then
 RETURN(FAIL):
fi:

if D1<>FAIL then
 D1:=D1-5:
fi:


lu:=min(seq(delt(L1[i]/L0[i],C1),i=90..100)):
lu:=evalf(lu,10):

[evalf(C1,30),identify(evalf(C1,30)),lu,D1]:


end:








#zeilL1rG(F,k,n,N,MaxC,r): Like zeilL1r(F,k,n,N,M,r), but also tries to guess a holonomic deserciption to the right side. Try:
#zeilL1rG(binomial(n,k)*binomial(n+k,k),k,n,N,40,2);
zeilL1rG:=proc(F,k,n,N,MaxC,r) local t,gu,i,L1,L2,ope,L0,M,opeR,D1:

option remember:

M:=max(seq((2+D1)*(1+MaxC-D1)+4,D1=0..MaxC)):
gu:=zeilL1(F,k,n,t,N,M):

if gu=FAIL then
  RETURN(FAIL):
fi:

ope:=gu[1]:


L0:=gu[2]:

L1:=gu[3]:

L2:=gu[4]:

L0:=subs(t=0,L1):

L1:= [seq(coeff(taylor(L1[i],t=0,r+6),t,r),i=1..nops(L1))]:
L2:= [seq(coeff(taylor(L2[i],t=0,r+6),t,r),i=1..nops(L2))]:

if L2=[0$nops(L2)] then 
RETURN( [ope,L0,L1,L2,[seq(evalf(L1[i]/L0[i],10),i=1..M)], [N-1,[0]]]):
fi:
opeR:=Findrec(L2,n,N,MaxC):

if opeR=FAIL then
 [ope,L0,L1,L2,[seq(evalf(L1[i]/L0[i],10),i=1..M)],FAIL]:
else
[ope,L0,L1,L2,[seq(evalf(L1[i]/L0[i],10),i=1..M)], [opeR,[op(1..degree(opeR,N),L2)]]]:
fi:
end:






#SeqFromRecG(ope,n,N,Ini,RHS1, K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)] and a holnomic description of the right side
#satisfied by the recurrence ope(n,N)f(n)=RHS1(n)
#extends it to the first K values. Try:
#SeqFromRecG(N-2,n,N,[1],[N-1,[4]],30);
SeqFromRecG:=proc(ope,n,N,Ini,RHS1,K) local yemin,coe1,i,L,ope1,gu,n1,j1,L1,L2:
yemin:=SeqFromRec(RHS1[1],n,N,RHS1[2],K):




coe1:=Yafe(ope,N)[1]:


yemin:=[seq(yemin[i]/subs(n=i,coe1),i=1..nops(yemin))]:


ope1:=Yafe(ope,N)[2]:


L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L) +yemin[n1-nops(Ini)] ]:
od:


L1:=ApplyOpe(ope,n,N,gu):
L2:=SeqFromRec(RHS1[1],n,N,RHS1[2],nops(L1)):

if not L1=L2 then
print(L1,L2):
 print(`Something went wrong`):
 RETURN(FAIL):
fi:
gu:
end:


#AZcG(F,k,n,r,MaxC,M):  Inputs a 
#hypergemetric term F (phrased in terms of the varlaibes k and n). It also inputs a small positive integer r, a medim-sized integer MaxC and
#a large integer M.
#It outputs the  approx. of the Apery-Zuidlin 
#constant of the coeff. of t^r the Zudilin-Straub transform of F, followed by the estimated delta. Try:
#AZcG(binomial(n,k)*binomial(n+k,k),k,n,2,15,1000);
AZcG:=proc(F,k,n,r,MaxC, M) local gu,N,ope,L0,L1,C1,C2,D1,Ini0,Ini1,i,lu,L1new:

if M<100 then
 print(M, `should have been at least 100`):
 RETURN(FAIL):
fi:

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


gu:=zeilL1rG(F,k,n,N,MaxC,r):




if gu=FAIL or gu[6]=FAIL then
 RETURN(FAIL):
fi:

if convert(gu[4],set)={0} then
RETURN(AZc(F,k,n,r,M)):
fi:


ope:=gu[1]:

L0:=gu[2]:

L1:=gu[3]:

Ini0:=[op(1..degree(ope,N),L0)]:
Ini1:=[op(1..degree(ope,N),L1)]:

L0:=SeqFromRec(ope,n,N,Ini0,M):

L1new:=SeqFromRecG(ope,n,N,Ini1,gu[6],M):

if L1new=FAIL then
 lprint([ope,n,N,Ini1,gu[6],M]):
 RETURN(FAIL):
fi:


if [op(1..nops(L1),L1new)]<>L1 then
 print(`Something is wrong`):

 RETURN(FAIL):
fi:

L1:=L1new:


C1:=evalf(L1[M]/L0[M]):
C2:=evalf(L1[M-1]/L0[M-1]):
 
D1:=trunc(-log[10](abs(C1-C2))):


if D1<20 then
 RETURN(FAIL):
fi:

D1:=D1-5:

if delt(L1[90]/L0[90],C1)=FAIL then
 RETURN(evalf(C1,30)):
fi:


lu:=min(seq(delt(L1[i]/L0[i],C1),i=90..100)):
lu:=evalf(lu,10):
[evalf(C1,30),identify(evalf(C1,30)),lu,D1]:


end:



#FindMaxkOld(F,n,k): Inputs a hypergeometric term in n and k, outputs the number such that k=alpha*n makes F(n,k) as large as possible. Try:
#FindMaxkOld(binomial(n,k),n,k);
FindMaxkOld:=proc(F,n,k) local F1,gu,i,lu,mu,ku;
F1:=convert(F,factorial):
gu:=normal(simplify(subs(k=k+1,F1)/F1)):


lu:=[solve(gu=1,k)]:


for i from 1 to nops(lu) do
mu:=lu[i]:

mu:=evalf(subs(n=1000,mu)):

mu:=coeff(mu,I,1):
if abs(mu)<10^(-20) then
 ku:=lu[i]:
 ku:=coeff(asympt(ku,n,4),n):
 if evalf(ku)>0 and evalf(ku)<1 then
 RETURN(ku):
 fi:
fi:
od:
FAIL:
end:


#FindMax(F,n,k): Inputs a hypergeometric term in n and k, outputs the number such that k=alpha*n makes F(n,k) as large as possible. Try:
#FindMaxk(binomial(n,k),n,k);
FindMaxk:=proc(F,n,k) local F1,gu,i,lu,mu,ku;
F1:=convert(F,factorial):
gu:=normal(simplify(subs(k=k+1,F1)/F1)):


lu:=[solve(gu=1,k)]:


for i from 1 to nops(lu) do
mu:=lu[i]:

mu:=evalf(subs(n=1000,mu)):

mu:=coeff(mu,I,1):
if abs(mu)<10^(-20) then
 ku:=lu[i]:
 ku:=op(1,asympt(ku,n,4))/n:
 if type(evalf(ku),numeric) and evalf(ku)>0 and evalf(ku)<1 then
 RETURN(ku):
 fi:
fi:
od:
FAIL:
end:


#cnk(n,k,r,L): The coefficient of the Taylor expansion of t^r in (1/((1+t)*...*(1+t/k))*((1-t)*...*(1-t/(n-k)))^L. Try:
#[seq(cnk(10,k,2,3),k=0..10)];
cnk:=proc(n,k,r,L) local gu,i,t:

gu:=mul(1+t/i,i=1..k)*mul(1-t/i,i=1..n-k):
gu:=1/gu^L:
gu:=taylor(gu,t=0,r+1):
coeff(gu,t,r):
end:


#cnk2(n,k,L): The coefficient of the Taylor expansion of t^2 in (1/((1+t)*...*(1+t/k))*((1-t)*...*(1-t/(n-k)))^L. Using the closed-form formula
#[seq(cnk2(10,k,3),k=0..10)];
cnk2:=proc(n,k,L) local i:
L/2*(add(1/i^2,i=1..k)+add(1/i^2,i=1..n-k))+
(1/2)*L^2*(-add(1/i,i=1..k)+add(1/i,i=1..n-k))^2:
end:


#Paper2(K1,K2): Inputs two positive integers K1 (at least 3) and K2 and a large integer M (for computing the Apery limit accurately), and outputs a paper with all the Apery limits gotten from
#considering the coefficient of t^2 in the Zudilin-Straub transform of binomial(n,k)^L*a^k for L from 3 to K1 and a from 1 to K2.
#Try:
#Paper2(4,2,2000);

Paper2:=proc(K1,K2,M) local gu,L,a,n,k,i,co,ope,F,N,ka,X,i1,A,B,lu,INI0,INI1,C1,C2,ka1:

print(`Lots of Interesting Apery Limits that arise from the coeff. of t^2 in the Zudilin-Straub transform of the binomial coefficients sum`):
print(``):
print(Sum(binomial(n,k)^L*a^k,k=0..n)):
print(``):
print(`for integers L from 3 to`, K1, `and  integers a from  1 to`, K2):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let c(n,k;L) be the following discrete function where 0<=k<=n and L is a positive integer`):
print(``):
print(L/2*(Sum(1/i^2,i=1..k)+Sum(1/i^2,i=1..n-k))+L^2/2*(Sum(-1/i,i=1..k)+ Sum(1/i,i=1..n-k))^2):
print(``):

co:=0:

for L from 3 to K1 do
print(`L is`, L):
 for  a from 1 to K2 do
print(` a is `, a):

 co:=co+1:
 F:=binomial(n,k)^L*a^k:
 ope:=zeil(F,k,n,N)[1]:
  gu:=zeilL1r(F,k,n,N,degree(ope,N)+3,2) :
  INI0:=[op(1..degree(ope,N),gu[2])]:
  INI1:=[op(1..degree(ope,N),gu[3])]:
  ka:=FindMaxk(F,n,k):
print(``):
print(`---------------------------------------------`):
print(`The next theorem is inspired by the coefficient of t^2 in the Zudilin-Straub transform of`, Sum(F,k=0..n)):
print(``):
 print(`Theorem Number`, co):
  print(` Let C be the limit, as n goes to infinity of the above c(n,k;L) with L=`,L, `and k= the integer part of`,ka*n, `or in floats`, evalf(ka,10)*n ):
    print(`then a  much faster way to compute this number is as the following Apery limit`):
 print(``):
 print(`Let A(n) and B(n) be the solution of the following linear recurrence `):
 print(``):
 print(add(coeff(ope,N,i1)*X(n+i1),i1=0..degree(ope,N))=0):
 print(``):
 print(`or in Maple format`):
 print(``):
 lprint(add(coeff(ope,N,i1)*X(n+i1),i1=0..degree(ope,N))=0):
 print(``):
 print(`with initial conditions`):
 print(``):
print( seq(A(i1)=INI1[i1],i1=1..nops(INI1))):
 print(``):
print( seq(B(i1)=INI0[i1],i1=1..nops(INI0))):
 print(``):
 print(`Then the sequence of quotients A(n)/B(n) converges very fast to C`):
 print(``):
 lu:=AZc(F,k,n,2,M) :
 print(``):
#print(`The error is estimated to be <=`, CONSTANT/B(n)^(1+lu[3])):

 print(``):
 print(` Using `, M, ` terms yield `, lu[4], `(decimal)  digits`):
 print(``):
 print(`Here it is to 30 digits`, lu[1]):
 print(``):
 print(`To illustrate how fast the convergence is, and also as check, let's compute the n=5000 and n=10000 values of the above sequence, whose limit is our C`):

ka1:=evalf(ka):
 C1:=evalf(
L/2*(add(1/i^2,i=1..trunc(ka1*5000))+add(1/i^2,i=1..5000-trunc(ka1*5000)))+ 
(L^2/2)*(-add(1/i,i=1..trunc(ka1*5000))+add(1/i,i=1..5000-trunc(ka1*5000)))^2,
30):


 C2:=evalf(
L/2*(add(1/i^2,i=1..trunc(ka1*10000))+add(1/i^2,i=1..10000-trunc(ka1*10000)))+ 
(L^2/2)*(-add(1/i,i=1..trunc(ka1*10000))+add(1/i,i=1..10000-trunc(ka1*10000)))^2,
30):

 print(``):
  print(C1,C2):
print(`as you can see the convergence is extremely slow using the definition`):
 print(``):
od:
od:

end:

#EvalAjnk(j,n,k): Evaluates A_j(n,k):=1/j*((-1)^(j-1)*Sum(1/i^j,i=1..k) - Sum(1/i^j,i=1..n-k) ). Try:
#EvalAjnk(5,20,5);
EvalAjnk:=proc(j,n,k) local i:

1/j*((-1)^(j-1)*add(1/i^j,i=1..k) - add(1/i^j,i=1..n-k) ):


end:

#CrL(A,r,L): The symbolic expression for  cnk(n,k,r,L) (q.v.) in terms of A[j]'s where it stands for A[j](n,k). Try:
#CrL(A,4,L);
CrL:=proc(A,r,L) local gu,j,t:
option remember:
gu:=add(A[j]*t^j,j=1..r):
gu:=exp(-L*gu):
gu:=taylor(gu,t=0,r+1):
expand(coeff(gu,t,r)):
end:




#cnkD(n,k,r,L): Same input and output as cnk(n,k,r,L) (q.v.) but via the explicit expression CrL(A,r,L) (q.v.). Try
#[seq(cnkD(10,k,2,3),k=0..10)];

cnkD:=proc(n,k,r,L) local gu,A,j:

gu:=CrL(A,r,L):

subs(seq(A[j]=EvalAjnk(j,n,k),j=1..r),gu):

end:




#EvalAjnk(j,n,k): Evaluates A_j(n,k):=1/j*((-1)^(j-1)*Sum(1/i^j,i=1..k) - Sum(1/i^j,i=1..n-k) ). Try:
#EvalAjnk(5,20,5);
EvalAjnk:=proc(j,n,k) local i:

1/j*((-1)^(j-1)*add(1/i^j,i=1..k) - add(1/i^j,i=1..n-k) ):


end:

#CrL(A,r,L): The symbolic expression for  cnk(n,k,r,L) (q.v.) in terms of A[j]'s where it stands for A[j](n,k). Try:
#CrL(A,4,L);
CrL:=proc(A,r,L) local gu,j,t:
option remember:
gu:=add(A[j]*t^j,j=1..r):
gu:=exp(-L*gu):
gu:=taylor(gu,t=0,r+1):
expand(coeff(gu,t,r)):
end:




#Paper(K1,K2): Inputs two positive integers K1 (at least 3) and K2 and a large integer M (for computing the Apery limit accurately), and outputs a paper with all the Apery limits gotten from
#considering the coefficient of t^r in the Zudilin-Straub transform of binomial(n,k)^L*a^k for L from 3 to K1 and a from 1 to K2 and r from 1 to L-1
#Try:
#Paper(4,2,2000);

Paper:=proc(K1,K2,M) local gu,L,a,n,k,i,co,ope,F,N,ka,X,i1,A,B,lu,INI0,INI1,C1,C2,ku,r,K,j,t:

print(`Lots of Interesting Apery Limits that arise from the coeff. of t^r in the Zudilin-Straub transform of the binomial coefficients sum`):
print(``):
print(Sum(binomial(n,k)^L*a^k,k=0..n)):
print(``):
print(`for integers L from 3 to`, K1, `and  integers a from  1 to`, K2, `and integers r from 2 to L-1`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`For positive integers k such that 0<=k<=n, and positive integer j, Let     `):
print(``):
print(K[j](n,k)=-1/j*((-1)^(j-1)*Sum(1/i^j,i=1..k)-Sum(1/i^j,i=1..n-k))):
print(``):

co:=0:

for L from 3 to K1 do
 for r from 2 to L do
 for  a from 1 to K2 do

if not (a=1  and r mod 2=1) then
 co:=co+1:
 F:=binomial(n,k)^L*a^k:
 lu:=AZc(F,k,n,r,M) :
if lu<>FAIL then
  ope:=zeil(F,k,n,N)[1]:
  gu:=zeilL1r(F,k,n,N,degree(ope,N)+3,r) :
  INI0:=[op(1..degree(ope,N),gu[2])]:
  INI1:=[op(1..degree(ope,N),gu[3])]:
  ka:=FindMaxk(F,n,k):
print(``):
print(`---------------------------------------------`):
print(``):
print(`The following theorem is inspired by  the coefficient of `, t^r , ` in the Zudilin-Straub t-transform of`):
print(``):
print(Sum(F,k=0..n)):
print(``):
print(``):
 print(`Theorem Number`, co):

  print(` Let C be the limit, as n goes to infinity of`):

  ku:=CrL(K,r,L):
 ku:=subs({seq(K[j]=K[j](n,k),j=1..r)},ku):

 print(ku):

print(`and k= the integer part of`,ka*n, `or in floats`, evalf(ka,10)*n ):
    print(`then a  much faster way to compute this number is as the following Apery limit`):
 print(``):
 print(`Let A(n) and B(n) be the solution of the following linear recurrence `):
 print(``):
 print(add(coeff(ope,N,i1)*X(n+i1),i1=0..degree(ope,N))=0):
 print(``):
 print(`or in Maple format`):
 print(``):
 lprint(add(coeff(ope,N,i1)*X(n+i1),i1=0..degree(ope,N))=0):
 print(``):
 print(`with initial conditions`):
 print(``):
print( seq(A(i1)=INI1[i1],i1=1..nops(INI1))):
 print(``):
print( seq(B(i1)=INI0[i1],i1=1..nops(INI0))):
 print(``):
 print(`Then the sequence of quotients A(n)/B(n) converges very fast to C`):
 print(``):
 lu:=AZc(F,k,n,r,M) :
 print(``):
 print(` Using `, M, ` terms yield `, lu[4], `(decimal)  digits`):
 print(``):
 print(`Here it is to 30 digits`, lu[1]):
 print(``):

 print(`To illustrate how fast the convergence is, and also as check, let's compute the n=5000 and n=10000 values of the above sequence, whose limit is our C`):

C1:=evalf(cnkD(5000,trunc(5000*evalf(ka)),r,L),20):
C2:=evalf(cnkD(10000,trunc(10000*evalf(ka)),r,L),20):
 print(``):
  print(C1,C2):
print(`as you can see the convergence is extremely slow using the definition`):
 print(``):

fi:
fi:
od:
od:
od:

end:


####new stuff


#cnkF(F,n,k,n1,k1,r): The coefficient of the Taylor expansion of t^r in  
#ZT(F,k,t,RF)/F for numeric n=n1 and k=k1. Try:
#cnkF(binomial(n,k)^3,n,k,5,2);

cnkF:=proc(F,n,k,n1,k1,r) local gu,t,RF:
gu:=ZT(F,k,t,RF)/F:
gu:=expand(normal(eval(subs({n=n1,k=k1},subs(RF=rf,gu))))):
gu:=taylor(gu,t=0,r+1):
coeff(gu,t,r):
end:


#EvalFrb(r,b,N): the numeric value of the atomic quantity F_{r,b}(N):=(-1)^r/r*Sum(b^r/j^r,j=1..N). Try:
#EvalFrb(2,3,10);
EvalFrb:=proc(r,b,N) local j:
(-1)^r/r*add(b^r/j^r,j=1..N):
end:



#EvalZTF(F,k,n,k1,n1,r): Evaluates the coeff. of t^r in the Taylor expansion of ZT(F,k,n,RF)/F at k=1,n=n1. Try
#EvalZTF(binomial(n,k)^3,k,n,10,20);
EvalZTF:=proc(F,k,n,k1,n1,r) local gu,RF,t:
gu:=normal(subs(RF(1,n)=n!,ZT(F,k,t,RF)/convert(F,factorial))):
gu:=subs(RF=rf,gu):
gu:=normal(simplify(eval(subs({n=n1,k=k1},gu)))):
coeff(taylor(gu,t=0,r+3),t,r):
end:


#PaperG(F,k,n,r,MaxC,M):Inputs a hypergeometric term F in discrete variables n and k, a small positive integer r, a medium integer MexC and a large M
#gives information about the Apery constant inspired by the coeff. of t^t in the Straub-Zudilin transform  if Sum(F,k=0..n). Try
#PaperG(binomial(n+k,k)*binomial(n,k),k,n,2,15,2000);
PaperG:=proc(F,k,n,r,MaxC,M) local gu1,gu2,t,N,ka,RF,yemin,A,B,C,c,R,ope,INI,opeY,IniY,C1,C2,t0,i:
t0:=time():
gu1:=zeilL1rG(F,k,n,N,MaxC,r):
if gu1=FAIL then
 RETURN(FAIL):
fi:

gu2:=AZcG(F,k,n,r,MaxC, M): 



if gu2=FAIL then
print(`MaxC=`, MaxC, `did not work out, try to raise it`):
 RETURN(FAIL):
fi:

  ka:=FindMaxk(F,n,k):

print(`The Apery Constant Inspired by the coefficient of `, t^r, ` in the Zudilin-Straub transform of`, Sum(F,k=0..n)):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let RF(t,k), as usual, be the raising factorial, t*(t+1)...(t+k-1) `):
print(``):
print(`We are interested in the limit, as n goes to infinity, of the coefficient of `, t^r, ` in the following quantity`):
print(``):
print(normal(subs(RF(1,n)=n!,ZT(F,k,t,RF)/convert(F,factorial)))):
print(``):
print(`and k equals`, ka*n, `that in floating-point is`, evalf(ka,20)*n):
print(``):
print(`Let's call this constant c (it can also be described directly in terms of harmonic-type expressions). Using this definition the convergence is extremely slow.`):
print(``):
print(`We will show how to compute it much faster, with exponential error-rate `):
print(``):

ope:=gu1[1]:
INI:=gu1[2]:
yemin:=gu1[6]:

print(`Let A(n) be the sequence of integers satisfying the linear recurrence`):
print(``):
print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Subject to the initial condition `):
print(``):
lprint(seq(A(i)=INI[i],i=1..degree(ope,N))):
print(``):

if yemin[2]=[0$nops(yemin[2])] then
print(``):
print(`Let  B(n) be the sequence satisfying the same recurrence i.e. `):
print(``):
lprint(add(coeff(ope,N,i)*B(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but with the initial conditions `):
print(``):
lprint(seq(B(i)=gu1[3][i],i=1..degree(ope,N))):

else

opeY:=yemin[1]:
IniY:=yemin[2]:

print(``):
print(`Let  B(n) be the sequence satisfying the INHOMOGENEOUS recurrence i.e. `):
print(``):
lprint(add(coeff(ope,N,i)*B(n+i),i=0..degree(ope,N))=C(n)):
print(``):
print(` with the initial conditions `):
print(``):
lprint(seq(B(i)=gu1[3][i],i=1..degree(ope,N)  )):
print(``):
print(`where the right side, C(n), satisfies the recurrence `):
print(``):
print(``):
lprint(add(coeff(opeY,N,i)*C(n+i),i=0..degree(opeY,N))=0):
print(``):
print(``):
print(` with the initial conditions `):
print(``):
lprint(seq(C(i)=IniY[i],i=1..nops(IniY))):
print(``):
print(`Note that it is very fast to compute many terms of C(n) and hence of B(n) `):
print(``):
fi:


 print(``):
#print(`The error is estimated to be <=`, CONSTANT/B(n)^(1+gu2[3])):
print(`Computing the`, M, `-th term gives `, gu2[4], ` digits `):


 print(``):
print(`The ratios B(n)/A(n) tend very fast to the constant c`):
print(``):
 print(`Here it is to 30 digits`, gu2[1]):
print(``):
if gu2[1]<>gu2[2] then
 print(`This is most probably`, gu2[2]):
print(``):
fi:


 print(`To illustrate how fast the convergence is, and also as check, let's compute the n=5000 and n=10000 values of the above sequence, whose limit is our c`):

C1:=evalf(EvalZTF(F,k,n,trunc(5000*evalf(ka)),5000,r),10):
C2:=evalf(EvalZTF(F,k,n,trunc(10000*evalf(ka)),10000,r),10):

 print(``):
  print(C1,C2):
print(`as you can see the convergence is extremely slow using the definition`):
 print(``):
print(`-------------------------------------------------------`):
 print(``):
print(`This ends this paper that took`, time()-t0, `seconds to generate. `):
 print(``):
end:








#PaperGshort(F,k,n,r,MaxC,M):Inputs a hypergeometric term F in discrete variables n and k, a small positive integer r, a medium integer MexC and a large M
#gives information about the Apery constant inspired by the coeff. of t^t in the Straub-Zudilin transform  if Sum(F,k=0..n). Try
#PaperG(binomial(n+k,k)*binomial(n,k),k,n,2,15,2000);
PaperGshort:=proc(F,k,n,r,MaxC,M) local gu1,gu2,t,N,ka,RF,yemin,A,B,C,c,R,ope,INI,opeY,IniY,C1,C2,t0,i:
t0:=time():
gu1:=zeilL1rG(F,k,n,N,MaxC,r):
if gu1=FAIL then
 RETURN(FAIL):
fi:

gu2:=AZcG(F,k,n,r,MaxC, M): 



if gu2=FAIL then
print(`MaxC=`, MaxC, `did not work out, try to raise it`):
 RETURN(FAIL):
fi:

  ka:=FindMaxk(F,n,k):

print(`The Apery Constant Inspired by the coefficient of `, t^r, ` in the Zudilin-Straub transform of`, Sum(F,k=0..n)):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let RF(t,k), as usual, be the raising factorial, t*(t+1)...(t+k-1) `):
print(``):
print(`We are interested in the limit, as n goes to infinity, of the coefficient of `, t^r, ` in the following quantity`):
print(``):
print(normal(subs(RF(1,n)=n!,ZT(F,k,t,RF)/convert(F,factorial)))):
print(``):
print(`and k equals`, ka*n, `that in floating-point is`, evalf(ka,20)*n):
print(``):
print(`Let's call this constant c (it can also be described directly in terms of harmonic-type expressions). Using this definition the convergence is extremely slow.`):
print(``):
print(`We will show how to compute it much faster, with exponential error-rate `):
print(``):

ope:=gu1[1]:
INI:=gu1[2]:
yemin:=gu1[6]:

print(`Let A(n) be the sequence of integers satisfying the linear recurrence`):
print(``):
print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0):
print(``):
print(`Subject to the initial condition `):
print(``):
lprint(seq(A(i)=INI[i],i=1..degree(ope,N))):
print(``):

if yemin[2]=[0$nops(yemin[2])] then
print(``):
print(`Let  B(n) be the sequence satisfying the same recurrence i.e. `):
print(``):
lprint(add(coeff(ope,N,i)*B(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but with the initial conditions `):
print(``):
lprint(seq(B(i)=gu1[3][i],i=1..degree(ope,N))):

else

opeY:=yemin[1]:
IniY:=yemin[2]:

print(``):
print(`Let  B(n) be the sequence satisfying the INHOMOGENEOUS recurrence i.e. `):
print(``):
lprint(add(coeff(ope,N,i)*B(n+i),i=0..degree(ope,N))=C(n)):
print(``):
print(` with the initial conditions `):
print(``):
lprint(seq(B(i)=gu1[3][i],i=1..degree(ope,N)  )):
print(``):
print(`where the right side, C(n), satisfies the recurrence `):
print(``):
print(``):
lprint(add(coeff(opeY,N,i)*C(n+i),i=0..degree(opeY,N))=0):
print(``):
print(``):
print(` with the initial conditions `):
print(``):
lprint(seq(C(i)=IniY[i],i=1..nops(IniY))):
print(``):
print(`Note that it is very fast to compute many terms of C(n) and hence of B(n) `):
print(``):
fi:


 print(``):
#print(`The error is estimated to be <=`, CONSTANT/B(n)^(1+gu2[3])):
print(`Computing the`, M, `-th term gives `, gu2[4], ` digits `):


 print(``):
print(`The ratios B(n)/A(n) tend very fast to the constant c`):
print(``):
 print(`Here it is to 30 digits`, gu2[1]):
print(``):
if gu2[1]<>gu2[2] then
 print(`This is most probably`, gu2[2]):
print(``):
fi:

 print(``):
print(`This ends this paper that took`, time()-t0, `seconds to generate. `):
 print(``):
end: