######################################################################
##BEUKERS.txt: Save this file as  BEUKERS.txt                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  BEUKERS.txt<Enter>                           #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Oct. 2019
 
print(`Created:  Oct. 2019`):
print(` This is BEUKERS.txt `):
print(`It is one of the Maple packages that accompany the article `):
print(` Automatic Discovery of Irrationailty Proofs and Irrationality Measures`):
print(`by Doron Zeilberger and Wadim Zudiln`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):



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

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

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

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):



ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: dn, Z2, Z2n, Z2F `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: deltEg, deltT, deltTS, Ope, Paper, Z2g, Z2gDelE, Z2gF, Z2gFS `):
 print(` `):


elif nops([args])=1 and op(1,[args])=deltEg then
print(`deltEg(f,X1,X2,C1,C2): Given f a linear combination of 1, X1 standing for the constant C1, and X2 standing for the constant C2,`):
print(` with integer coefficients, finds`):
print(`the empirical delta such that abs(f)=1/max(coeff(f,1),coeff(f,X1),coeff(f,X2))^delta. Try:`):
print(`deltEg(numer(Z2gFS(20,3,X1,X2)),X1,X2,dilog(2/3), log(2/3));`):

elif nops([args])=1 and op(1,[args])=deltT then
print(`deltT(a): The theoretical delta such that |A(n)+B(n)dilog((a-1)/a)+C(n)log(a/(a-1))|<=C/max(A(n),B(n),C(n))^delta.`):
print(`where A(n),B(n),C(n) are the integer sequences gotten from the generalized Beukers integral`):
print(`int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1), x=0..1),y=0..1)=A1(n)+B1(n)*dilog(a/(a-1))+C1(n)*log(a/(a-1))`):
print(`where A1(n), B1(n), C1(n) are RATIONAL numbers, but defining`):
print(`A(n)=lcm(1..n)^2*A1(n), B(n)=lcm(1..n)^2*B1(n), C(n)=lcm(1..n)^2*C1(n), gives you integer sequences . Try:`):
print(`deltT(3);`):


elif nops([args])=1 and op(1,[args])=deltTS then
print(`deltTS(a): Like deltT(a) but for SYMBOLIC a. Try:`):
print(`deltTS(a);`):


elif nops([args])=1 and op(1,[args])=dn then
print(`dn(n): lcm(1...n). Try: dn(10);`):


elif nops([args])=1 and op(1,[args])=Ope then
print(`Ope(n,N,a): The recurrence operator in n, where N is the forward shift operator, annihilating  int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1). `):
print(`It was obtained using the Apagodu-Zeilberger multivariable extension of the ALmkvist-Zeilberger algorithm.`):
print(`Try: `):
print(`Ope(n,N,a);`):

elif nops([args])=1 and op(1,[args])=Paper then
print(`Paper(K):  A computer generated paper about the joint approximation of DiLog((a-1)/a) and Log((a-1)/a) for ANY integer a>=2. Try:`):
print(`K is a positve integer telling how many examples of delta to give. Try:`):
print(`Paper(5): `):

elif nops([args])=1 and op(1,[args])=Z2 then
print(`Z2(n): The Beukers integral used in his proof that Zeta(2) is irrational. Done Directly. Very slow. You should use Z2F(n,X). try:`):
print(`[seq(Z2(i),i=1..20)];`):



elif nops([args])=1 and op(1,[args])=Z2F then
print(`Z2F(n,X): The Beukers integral used in his proof that Zeta(2) is irrational. Here X is a symbol that stands for Zeta(2)=Pi^2/6`):
print(`It is done Via the second-order recurrence (gotten from the Almkvist-Zeilbeger algorithm). MUCH FASTER. try:`):
print(`[seq(Z2F(i,Pi^2/6),i=1..100)];`):


elif nops([args])=1 and op(1,[args])=Z2n then
print(`Z2n(n): The n-th approximation to Zeta(2) implied by the Beukers integral. Try:`):
print(`[seq(Z2n(n),n=1..100)]`):

elif nops([args])=1 and op(1,[args])=Z2g then
print(`Z2g(n,a): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1),x=0..1),y=0..1);`):
print(`Done DIRECTLY. Very slow. For fast computation, use Z2gF(n,a). Try`):
print(`[seq(Z2g(i,3),i=1..10)];`):


elif nops([args])=1 and op(1,[args])=Z2gDelE then
print(`Z2gDelE(a,N): The smallest empirical delta for the linear forms in dilog(a/(a-1)) and log(a/(a-1)) from N to 2N. Try:`):
print(`Z2gDelE(3,100);`):


elif nops([args])=1 and op(1,[args])=Z2gF then
print(`Z2gF(n,a): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1),x=0..1),y=0..1);`):
print(`Done via the third-order linear recurrence gotten from the Apagodu-Zeilberger multivariable extension of theAlmkvist-Zeilberger algorithm .Very Fast. Try`):
print(`[seq(Z2gF(i,3),i=1..100)];`):


elif nops([args])=1 and op(1,[args])=Z2gFS then
print(`Z2gFS(n,a,X1,X2): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1),x=0..1),y=0..1);`):
print(`BUT now X1 and X2 are SYMBOLS that stand for dilog((a-1)/a) and log((a-1)/a) . It is for convenience for extraction the thee sequences.`):
print(`It is done via the third-order linear recurrence gotten from the Apagodu-Zeilberger multivariable extension of theAlmkvist-Zeilberger algorithm .Very Fast. Try`):
print(`[seq(Z2gFS(i,3,X1,X2),i=1..100)];`):

else
print(`There is no ezra for`,args):
fi:
 
end:

Digits:=7000:

#########From EKHAD

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

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

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


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

AZd1:= proc(A,y,n,N) local gu,ORDER,i1:
option remember:
gu:=AZd(A,y,n,N):

if gu<>0 then
gu:=gu[1]:
ORDER:=degree(gu,N):
gu:=gu/coeff(gu,N,ORDER):
gu:=1-subs(n=n-ORDER,gu)/N^ORDER:
gu:=add(factor(coeff(gu,N,-i1))*N^(-i1),i1=1..ORDER):
 RETURN(gu):
fi:

0:
end:

#########end From EKHAD


#dn(n): lcm(1...n). Try: dn(10);
dn:=proc(n) local i:
lcm(seq(i,i=1..n)):
end:


ez:=proc(): print(`     deltEg(f,X1,X2,C1,C2),  `): 
print(`Z2gDelE(a,N)`):
print(`deltEg(numer(Z2gFS(20,3,X1,X2)),X1,X2,dilog(2/3), log(2/3));, Ope(n,N,a), Alpha(a), deltT(a)  `):

end:


dn:=proc(n) local i:lcm(seq(i,i=1..n)):end:

Digits:=3000:

#deltE(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)):
deltE:=proc(a,c): 
if type(a,integer) then
 RETURN(FAIL):
fi:

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

#Z2F(n,X) int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y)^(n+1),x=0..1),y=0..1), where X=zeta(2) (aka =Pi^2/6)
Z2F:=proc(n,X) option remember:

if n=0 then
 X:
elif n=1 then
 5-3*X:
else
-(11*n^2-11*n+3)/n^2*Z2F(n-1,X)+(-1+n)^2/n^2*Z2F(n-2,X):
fi:
end:


#Z2(n): The Beukers integral, try:
#Z2(10);
Z2:=proc(n) local F,x,y:
F:=x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y)^(n+1):
int(int(F,x=0..1),y=0..1):
end:


##This part, about Zeta 3, does not work in Maple hence commented
#Z3(n): The Beukers integral for Zeta(3). Done directly. Try:
#Z3(5);
#Z3:=proc(n) local F,x,y,w:
#F:=x^n*(1-x)^n*y^n*(1-y)^n*w^n*(1-w)^n/(1-(1-x*y)*w)^(n+1):
#int(int(int(F,x=0..1),y=0..1),w=0..1):
#end:

#Z3fl(n):The Beukers integral for Zeta(3). Done directly. Try: Floating point
#Z3fl:=proc(n) local F,x,y,w:
#F:=x^n*(1-x)^n*y^n*(1-y)^n*w^n*(1-w)^n/(1-(1-x*y)*w)^(n+1):
#evalf(Int(Int(Int(F,x=0..y),y=0..w),w=0..1)):
#end:
#End of part about Zeta 3, does not work in Maple hence commented

#Z2n(n): The n-th approximation of zeta(2)
Z2n:=proc(n) local gu,X:
gu:=Z2F(n,X):
-coeff(gu,X,0)/coeff(gu,X,1):
end:



##Z2G2(n,a,b):int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x^a*y^b)^(n+1),x=0..1),y=0..1):
#Z2G2:=proc(n,a,b) local F,x,y:
#F:=x^n*(1-x)^n*y^n*(1-y)^n/(1-x^a*y^b)^(n+1):
#int(int(F,x=0..1),y=0..1):
#end:




#Z2g(n,a): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1),x=0..1),y=0..1);
Z2g:=proc(n,a) local F,x,y:
F:=x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1):
int(int(F,x=0..1),y=0..1):
end:



#Z2gF(n,a): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1),x=0..1),y=0..1) The fast way, using the recurrence
Z2gF:=proc(n,a) option remember:

if n=0 then
 a*dilog((a-1)/a):
elif n=1 then
-a^2*(2*a+1)*dilog((a-1)/a)+3*(a-1)*a^2*log((a-1)/a)+5*a^2;
elif n=2 then
a^3*(6*a^2+12*a+1)*dilog((a-1)/a)-a^3*(29*a+9)*(a-1)/2*log((a-1)/a)-a^3*(82*a+43)/4:
else

-a*(256*a^2*n^3-672*a^2*n^2-120*a*n^3+496*a^2*n+300*a*n^2-81*n^3-120*a^2-230*a*n+207*n^2+60*a-141*n+30)/n^2/(32*a*n-52*a-27*n+42)*Z2gF(n-1,a)
-a^2*(512*a^3*n^3-1856*a^3*n^2-1072*a^2*n^3+2112*a^3*n+3856*a^2*n^2+636*a*n^3-720*a^3-4332*a^2*n-2268*a*n^2-81*n^3+1440*a^2+2504*a*n+288*n^2-800*a-309*n+90)/n^2/(32*a*n-52*a-27*n+42)*Z2gF(n-2,a)
+(-2+n)^2*(a-1)*(32*a*n-20*a-27*n+15)*a^3/n^2/(32*a*n-52*a-27*n+42)*Z2gF(n-3,a):

fi:

end:


#Z2gFS(n,a,X1,X2): Like Z2gF(n,a) but in terms of X1 and X2 stading for dilog(a/(a-1)) and log(a/(a-1)) respectively. Try:
#Z2gFS(10,3,X1,X2);
Z2gFS:=proc(n,a,X1,X2) option remember:

if n=0 then
 a*X1:
elif n=1 then
-a^2*(2*a+1)*X1+3*(a-1)*a^2*X2+5*a^2;
elif n=2 then
a^3*(6*a^2+12*a+1)*X1-a^3*(29*a+9)*(a-1)/2*X2-a^3*(82*a+43)/4:
else

-a*(256*a^2*n^3-672*a^2*n^2-120*a*n^3+496*a^2*n+300*a*n^2-81*n^3-120*a^2-230*a*n+207*n^2+60*a-141*n+30)/n^2/(32*a*n-52*a-27*n+42)*Z2gFS(n-1,a,X1,X2)
-a^2*(512*a^3*n^3-1856*a^3*n^2-1072*a^2*n^3+2112*a^3*n+3856*a^2*n^2+636*a*n^3-720*a^3-4332*a^2*n-2268*a*n^2-81*n^3+1440*a^2+2504*a*n+288*n^2-800*a-309*n+90)/n^2/(32*a*n-52*a-27*n+42)*Z2gFS(n-2,a,X1,X2)
+(-2+n)^2*(a-1)*(32*a*n-20*a-27*n+15)*a^3/n^2/(32*a*n-52*a-27*n+42)*Z2gFS(n-3,a,X1,X2):

fi:

end:









#deltEg(f,X1,X2,C1,C2): Given f a linear combination of 1, X1 standing for the constant C1, and X2 standing for the constant C2,
# with integer coefficients, finds
#the empirical delta such that abs(f)=1/max(coeff(f,1),coeff(f,X1),coeff(f,X2))^delta. Try:
#deltEg(numer(Z2gFS(20,3,X1,X2)),X1,X2,dilog(2/3), log(2/3));
deltEg:=proc(f,X1,X2,C1,C2) local f1,q,p1,p2,L,E1:

f1:=numer(f):
q:=coeff(coeff(f1,X1,0),X2,0):
p1:=coeff(f1,X1,1):
p2:=coeff(f1,X2,1):
L:=max(abs(q),abs(p1),abs(p2)):

E1:=evalf(subs({X1=C1,X2=C2},f1)):


evalf(-log(abs(E1))/log(L)):

end:




#Z2gDelE(a,N): The smallest empirical delta for the linear forms in dilog(a/(a-1)) and log(a/(a-1)) from N to 2N. Try:
#Z2gDelE(3,100);
Z2gDelE:=proc(a,N) local X1,X2,n,gu:

gu:=min(seq(deltEg(numer(Z2gFS(n,a,X1,X2)),X1,X2,dilog((a-1)/a), log((a-1)/a)),n=N..2*N)):
evalf(gu,20):
end:



#Ope(n,N,a): The recurrence operator annihilating  int(x^n*(1-x)^n*y^n*(1-y)^n/(1-x*y/a)^(n+1). Try:
#Ope(n,N,a);
Ope:=proc(n,N,a)
-a^3*(1+n)^2*(a-1)*(32*a*n+76*a-27*n-66)/(32*a*n+44*a-27*n-39)/(3+n)^2+a^2*(512*a^3*n^3+2752*a^3*n^2-1072*a^2*n^3+4800*a^3*n-5792*a^2*n^2+636*a*n^3+2736*a^3-10140*a
^2*n+3456*a*n^2-81*n^3-5796*a^2+6068*a*n-441*n^2+3472*a-768*n-432)/(32*a*n+44*a-27*n-39)/(3+n)^2*N+a*(256*a^2*n^3+1632*a^2*n^2-120*a*n^3+3376*a^2*n-780*a*n^2-81*n^3
+2232*a^2-1670*a*n-522*n^2-1170*a-1086*n-717)/(32*a*n+44*a-27*n-39)/(3+n)^2*N^2+N^3:
end:

Ope1:=proc(N,a) local n,ope,i:
ope:=op(2,factor(lcoeff(numer(Ope(n,N,a)),n))):
add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
end:


#Alpha(a): The rate of growth of the sequences A(n),B(n),C(n) and the reciprocal of E(n,a). 
Alpha:=proc(a) local ope,n,N:
ope:=Ope(n,N,a):
ope:=numer(ope):
ope:=lcoeff(ope,n):
[solve(ope,N)]:
end:






#deltT(a): The theoretical delta such that |A(n)+B(n)dilog((a-1)/a)+C(n)log(a/(a-1))|<=C/max(A(n),B(n),C(n))^delta. Try
#deltT(3);
deltT:=proc(a) local gu,N,ka,ga,i:
gu:=-a^3*(a-1)+a^2*(16*a^2-20*a+3)*N+a*(8*a+3)*N^2+N^3:
gu:=evalf([solve(gu,N)]):

gu:=[seq(abs(coeff(gu[i],I,0)),i=1..nops(gu))]:

ka:=min(op(gu)):
ga:=max(op(gu)):

[ka,ga]:

-(log(ka)+2)/(log(ga)+2):
end:



#deltTS(a): The theoretical delta such that |A(n)+B(n)dilog((a-1)/a)+C(n)log(a/(a-1))|<=C/max(A(n),B(n),C(n))^delta. 
#For symbolic a. Try:
#deltTS(a);
deltTS:=proc(a) local gu,N,ka,ga:

gu:=-a^3*(a-1)+a^2*(16*a^2-20*a+3)*N+a*(8*a+3)*N^2+N^3:
gu:=[solve(gu,N)]:

ka:=gu[1]:
ga:=-gu[2]:
-(log(ka)+2)/(log(ga)+2):
end:


#ZZ(a): The rate of growth of E(n,a) and of the sequence A(n), B(n), C(n) for dilog((a-1)/a) and log((a-1)/a)
#For symbolic a. Try:
#ZZ(a);
ZZ:=proc(a) local gu,N,ka,ga:

gu:=-a^3*(a-1)+a^2*(16*a^2-20*a+3)*N+a*(8*a+3)*N^2+N^3:
gu:=[solve(gu,N)]:

ka:=gu[1]:
ga:=-gu[2]:
[ka,ga]:
end:


#Paper(K): outputs a paper about joint diophantine approximations of dilog((a-1)/a) and log(a/(a-1))
#with a proof (except for one lemma left to the reader). Try:
#K is a positve integer telling how many examples of delta to give
#Paper(3);
Paper:=proc(K) local a,n,delta,gu,a1,lu,f,vu:
end:
Paper:=proc(K) local a,gu,a1,E,n,x,y,lu,f,ku,vu,N,i:
gu:=deltTS(a):


print(`Integer Linear Combinations of Dilog((a-1)/a) and Log((a-1)/a) that are VERY close to an INTEGER`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Fix a  positive integer larger than 1, let's call it a. There exist EXPLICIT integer sequences A(n), B(n), C(n)`):
print(`(that depend on a, of course) such that `);
print(``):
print(`abs(A(n)+B(n)*DiLog((a-1)/a)+C(n)*Log((a-1)/a))<= CONSTANT/max(A(n),B(n),C(n))^delta, where`):
print(`delta equals `):
print(``):
print(gu):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu):
print(``):
print(`Before going into the proof, for the sake of conreteness, let state the  values for a from 2 to 20.`):
print(``):
print(` We also state the smallest empirical delta between n=100 and n=200 for these sequences to be described shortly`):
print(``):
for a1 from 2 to 1+K do
print(`a= `, a1, `delta=`, evalf(coeff(evalf(subs(a=a1,gu)),I,0),20),  `Smallest Empirical: `, Z2gDelE(a1,100)) :
od:
lu:=ZZ(a):
print(``):
print(`Proof of the Theorem`):
print(``):
print(`Consider a generalization of Beukers's famous integral that established the improved Apery proof of the irrationality of Zeta(2).`):
print(``):
print(`Let's call it`, E(n,a) ):
print(``):
print(E(n,a)=Int(Int(x^n*(1-x^n)*y^n*(1-y)^n/((1-x*y/a))^(n+1),x=0..1),y=0..1)):
print(``):
print(`It is readily seen that E(n,a) is a linear combination with RATIONAL number coefficients of 1, Dilog((a-1)/a)) and Log((a-1)/a)`):
print(``):
print(`E(n,A)=A1(n)+B1(n)*DiLog((a-1)/a)+ C1(n)*Log((a-1)/a) `):
eprint(``):
print(`This introduces three sequences of RATIONAL numbers, A1(n), B1(n), C1(n) (we supress the dependance on a for notational ease)`):
print(``):
print(`Lemma 1: The absolute value of E(n,a)  is asymptotic (up to a constant) to`):
print(``):
print(lu[1]^n):
print(``):
print(`Proof: Maximize the function`, x*(1-x)*y*(1-y)/(1-x*y/a), `in 0<=x,y<=1 using  two-variable calculus`):
print(``):
print(`Indeed taking partial derivatives with respect to x and y, and setting it equal to zero gives that the maximum is attained at the point`):
f:= x*(1-x)*y*(1-y)/(1-x*y/a):
ku:=solve({numer(diff(f,x)),numer(diff(f,y))},{x,y})[7]:
print(ku):
print(``):
print(`Now plug it into`,  x*(1-x)*y*(1-y)/(1-x*y/a) ):
print(``):
print(`Lemma 2: The three sequences of rational numbers A1(n), B1(n), C1(n) all are asymptotic, up to  polynomial corrections to `):
print(``):
print(lu[2]^n):
print(``):
print(`Proof: Using the amazing Multivariable extension of the Almkvist-Zeilberger algorithm, due to Moa Apagodu and Doron Zeilberger,`):
print(`it was discovered and proved that E(n,a) satisfies the following THIRD-ORDER linear recurrence  equation with polynomial coefficients`):
vu:=Ope(n,N,a):
print(``):
print(add(coeff(vu,N,i)*E(n+i,a),i=0..degree(vu,N))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(add(coeff(vu,N,i)*E(n+i,a),i=0..degree(vu,N))=0):
print(``):
print(`It follows that the three sequences of rational numbers, A1(n), B1(n), C1(n) also satisy this recurrence, so let denote them`):
print(`all by X(n)`):
print(``):
print(add(coeff(vu,N,i)*X(n+i),i=0..degree(vu,N))=0):
print(``):

print(`where X(n) stand for each of A1(n),B1(n), C1(n), that of course have different initial conditions. `):
print(``):
print(`Taking the leading term in n, gives us the constant-coefficient recurrence that approximates (up to polynomial corrections)`):
print(``):
print(`The sequences A1(n), B1(n), C1(n) `):
print(``):
vu:=Ope1(N,a):
print(``):
lprint(add(coeff(vu,N,i)*X1(n+i),i=0..degree(vu,N))=0):
print(``):
print(`where X1(n) is a C-finite approximation with the same asymptotics (up to polynomial factors), thanks to the Poincare lemma.`):
print(``):
print(`Solving the indicial equation `):
print(``):
lprint(add(coeff(vu,N,i)*x^i,i=0..degree(vu,N))=0):
print(``):
print(`and taking the largest root, establishes the lemma.`):
print(``):
print(`Lemma 3: Let d(n) be the least common multiple of the first n natural numbers. Recall that d(n) is asymptotic to exp(n)`):
print(``):
print(`A(n)=A1(n)*d(n)^2, B(n)=B1(n)*d(n)^2, C(n)=C1(n)*d(n)^2 are INTEGER sequencs `):
print(``):
print(`Proof: Left to the reader`):
print(``):
print(`We are now ready to prove the theorem: A(n), B(n), C(n) are asymptotic to`):
print(``):
print(lu[2]^n*exp(2*n)):
print(``):
print(`d(n)^2*E(n,a) is asymptotic to`):
print(``):
print(lu[1]^n*exp(2*n)):
print(``):
print(`Our delta is such that `):
print(``):
print(`(Max(A1(n),B1(n),C1(n))*exp(2*n))^delta=exp(2*n)/E(n,a)`):
print(``):
print(`Taking logarithms of both sides, and solving for delta, establishes the theorem. QED.`):
end: