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

Digits:=20000:

print(`Created: Dec. 2019`):
print(`  RUKHADZE.txt `):
print(`It is  one of the Maple package that accompanies the article `):
print(`Automatic Discovery of Irrationality Proofs and Irrationality Measures`):
print(`by  Doron Zeilberger and Wadim Zudiln`):
print(`available from the arxiv and from `):
print(` http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/gat.html `):
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(`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:  AZd, BU, EmpDelABCa, GuessK, GuessK1, GuessKr, GuessKrPair, MakeDB, NorPair, Nu1Extra`):
 print(` OpeABCa, OpeABCa1, SeqFromRec, SeqPair, ShoresABCa `):
 print(`  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  BestABCa, BestABCaG, deltE, InfoABCa, LABCa, PaperABCa, SeqABCa  `):
 print(` `):




elif nops([args])=1 and op(1,[args])=AZd then
print(`AZd(A,y,n,N): The Almkvist-Zeilberger algorithm. Inputs a hyper-exponential expression A in the continuous variable`):
print(`y and the discrete variable n, and and a symbol N, outputs the recurrence operator in forward shift N `):
print(`The f(n):=Int( A(y,n),y); (assuming a contour integral or one that vanishes at the end points`):
print(`satisfies: ope(n,N)f(n)=0. Try:`):
print(`AZd(x^n*(1-x)^n,x,n,N);`):


elif nargs=1 and args[1]=BestABCa then
print(`BestABCa(a,K,M): Explores the empirical delta for the terms between n=M-10 and n=M, in the rational approximations in the Rukhazde integral`):
print(`LABCa(A,B,C,a,n1) for 1<=A,B,C<=K`):
print(`it only records those with delta >cu `):
print(`BestABCa(1,3,500);`):

elif nargs=1 and args[1]=BestABCaG then
print(`BestABCaG(a,K,M,G): Explores the empirical delta for the terms between n=M-10 and n=M, in the rational approximations in the Rukhazde integral`):
print(`LABCa(A,B,C,a,n1) for 1<=A,B,C<=K and A and B are beween C-G and C+G whose delta is above that for (1,1,1). Try: `):
print(`BestABCaG(1,3,300,1);`):

elif nargs=1 and args[1]=BU then
print(`BU(f,a): If f is A+B*log((a+1)/a)) with A and B rational functions returns [A,B]`):
print(`Try: `):
print(`BU(11/19*log(5/4)-13/109);`):

elif nargs=1 and args[1]=EmpDelABCa then
print(`EmpDelABCa(A,B,C,a,M): the empirical deltas from M-10 to M for the Diophatine approximations implied by`):
print(`LABCa(A,B,C,a,n0) for n0 from M-10 to M. Try:`):
print(`EmpDelABCa(6,8,7,1,100);`):

elif nargs=1 and args[1]=deltE then
print(`deltE(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]=GuessK then
print(`GuessK(L,m,r): inputs a sequence  of rational numbers, L, a positive integer m, and a positive integer r`):
print(` finds a rational number K such that L[i]*K^i*dn(r*i) are all integers. Try:`):
print(`GuessK([seq((2/3)^i/dn(5*i),i=1..20)],30,5);`):

elif nargs=1 and args[1]=GuessK1 then
print(`GuessK1(L,m): Inputs a sequence of integers L, `):
print(` finds a constant K of the form K=ithprime(i)^j, with i and j<=m   such that L[n]/K^n is an integer for all n. Try:`):
print(`GuessK1([seq((i^3+1)*25^(i),i=1..30)],10);`):

elif nargs=1 and args[1]=GuessKr then
print(`GuessKr(L,m): inputs a sequence  of rational numbers, L, a positive integer m,`):
print(`finds a positive integer r, and a  rational number K [with search space given by m, see GuessK1(L,m)], `):
print(`such that L[i]*K^i*dn(r*i) are all integers. Try:`):
print(`GuessKr([seq((2/3)^i/dn(5*i),i=1..20)],30);`):

elif nargs=1 and args[1]=GuessKrPair then
print(`GuessKrPair(P,m): Given a pair,P, of sequences of rational numbers and a positive integer m`):
print(`finds a non-negative integer r and a rational number K such that`):
print(`BOTH P[1][n]*dn(r*n)*K^(i*n) and P[1][n]*dn(r*n)*K^(i*n) are sequences of INTEGERS. `):
print(`It tries out all the primes up to m and all the powers up to m.`):
print(`Type: `):
print(`GuessKrPair(SeqPair(1,1,1,1,40),20);`):


elif nargs=1 and args[1]=InfoABCa then
print(`InfoABCa(A,B,C,a,x,n,N,M): inputs positive integers A, B, C, and a such that gcd(A,B,C), symbols x,n,N, and a large positive integer M`):
print(`outputs a list consisting of `):
print(`(i) The integrand `):
print(`(ii) the Empirical delta for the approximation to Pi from the sequence obtained from the Rukhazde integral`):
print(`with A, B, C, a`):
print(`(iii) the linear recurrence operator in n and the shift operator N annihilating E(n) and hence the coefficients of log((a+1)/a) and 1`):
print(`(iv) The leading coefficient in n, in other words the constant-coefficient operator approximinating it.`):
print(`the so-called indicial equation`):
print(`(v) The pair [larger root, absolute value of smaller root]`):
print(`(vi) The conjectured pair [r,K] such that lcm(1..r*n)*K^n*E(n) has integer coefficients`):
print(`In other words C1(n):=lcm(1..r*n)*K^n*C(n) and  D1(n):=lcm(1..r*n)*K^n*D(n)  are INTEGERS sequences`):
print(`(vii) The upper bound for the irrationality measure not taking advantage of the fact that gcd(C1(n),D1(n)) may grow. `):
print(`(viii) the smallest  log(gcd(A1(n),A2(n))/n for n between M/2 and M`):
print(`(ix) assuming that this is a lower bound for the lim sup of  log(gcd(C1(n),D1(n))/n as n goes to infinity ,`):
print(`the implied upper bound for the irrationality measure of Pi. Try:`):
print(`InfoABCa(6,8,7,1,x,n,N,400);`):



elif nargs=1 and args[1]=LABCa then
print(`LABCa(A,B,C,a,n1): int((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x=0..1). Try:`):
print(`[seq(LABCa(1,1,1,1,n1),n1=0..10)];`):
print(`[seq(LABCa(6,8,7,1,n1),n1=0..10)];`):

elif nargs=1 and args[1]=MakeDB then
print(`MakeDB(M,n,N): inputs a positive integer M, and outputs a database of all the operators for computing`):
print(`LABCa(A,B,C,a,n) for 1<=A,B,C<=M and gcd(A,B,C)=1. The output is a table such that T[A,B,C] is the desired operator.`):
print(` Try: `):
print(`MakeDB(2,n,N) ; `):


elif nargs=1 and args[1]=NorPair then
print(` NorPair(P,r,K): Given a pair of sequences of rational numbers P, multiplies the n-th term by `):
print(` by lcm(1..r*n)*K^n. It also returns the list of gcds, let's call it G, and the list of log(G[i])/i and `):
print(` and the smallest and largest log(G[i])/i in the second half, in floating point. Try: `):
print(`NorPair(SeqPair(1,1,1,1,40),1,1);`):
print(`NorPair(SeqPair(6,8,7,1,40),7,1);`):


elif nargs=1 and args[1]=Nu1 then
print(`Nu1(a,b,K,r): if you have a sequence of diophantine approximations A(n)+B(n)*c such that`):
print(`(i)|A(n)+B(n)*c|<=Cb^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers`):
print(`outputs the implied upper bound for the irrationality measure.`):
print(`Here we do not use the fact that gcd(A(n),B(n)) may be of exponential growth.`):
print(`For the more general case, try Nu1Extra(a,b,K,r,Extra); `):
print(`Nu1(op(ShoreshPiAB(N,3,5)[2]),25/16,10);`):
print(`Nu1(op(ShoreshPiAB(N,2,3)[2]), 1/32,8);`):

elif nargs=1 and args[1]=Nu1Extra then
print(`Nu1Extra(a,b,K,r,K1): if you have a sequence of diophantine approximations A(n)+B(n)*c such that`):
print(`(i)|A(n)+B(n)*c|<=CONST*b^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers`):
print(`AND it looks like, empirically or rigorously that gcd(K^n*dn(r*n)*A(n), K^n*dn(r*n)*B(n)) are O(e^(K1*n))`):
print(`outputs the implied upper bound for the irrationality measure. Try:`):
print(`Nu1Extra(op(ShoreshPiAB(N,3,5)[2]),25/16,10,0);`):
print(`Nu1Extra(op(ShoreshPiAB(N,2,3)[2]),1/32,8,0);`):


elif nargs=1 and args[1]=OpeABCa then
print(`OpeABCa(A,B,C,a,n,N): the linear difference operator in backward notation for computing LABCa(A,B,C,a,n1). Try`):
print(`OpeABCa(6,8,7,1,n,N); `):


elif nargs=1 and args[1]=OpeABCa1 then
print(`OpeABCa1(A,B,C,a,n,N): the linear difference operator in backward notation for computing LABCa(A,B,C,a,n1), but in BACKWARD notation. Try`):
print(`OpeABCa1(6,8,7,1,n,N); `):

elif nargs=1 and args[1]=PaperABCa then
print(`PaperABCa(A,B,C,a,M): inputs relatively prime positive integers A, B, C, and a positive integer a, and a large positive integer M, at least 100`):
print(`outputs an article about diophantine approximations to log((a+1)/a) implied`):
print(`the Rukhadze integral LABCa(A,B,C,a,n,x). `):
print(`as given in Wadim Zudlin's paper, An essay on irrationality measures of pi and other logarithms, arXiV: math/0404523v2 [math.NT] 27 May 2004`):
print(`bottom of page 4`):
print(`It proves rigorously (modulo a simple divisibility lemma) a coarse upper bound for the irrationality measure  of log((a+1)/a). Try:`):
print(`PaperABCa(1,1,1,1,400);`):
print(`PaperABCa(6,8,7,1,400);`):

elif nargs=1 and args[1]=SeqABCa then
print(`SeqABCa(A,B,C,a,M): The first M terms of  LABCa(A,B,C,a,n0) for n0 from 1 to M. `):
print(`SeqABCa(1,1,1,1,100);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,INI,M): given a linear recurrence operator ope in n and the negative shift operator N^(-1) and`):
print( ` (a list of initial values (starting at n=0), outputs the first  M+1 terms of the sequence starting at n=0. `):
pritn(` Try: `):
print(` SeqFromRec(1/N+1/N^2,n,N,[0,1],10); `):

elif nargs=1 and args[1]=SeqPair then
print(`SeqPair(A,B,C,a,M): outputs a pair of integer sequences of length M.`):
print(`The first one is the first M terms, starting at n=1, of the coefficient of 1 in`):
print(`LABCa(A,B,C,a,n1) and the second the coefficient of log(a+1). Try:`):
print(`SeqPair(6,8,7,1,100);`):

elif nargs=1 and args[1]=ShoreshABCa then
print(`ShoreshABCa(N,A,B,C,a): the indicial polynomial of the recurrence for the Rukhadze sequence for (A,B,C,a)`):
print(`absolute values of the largest and smallest roots of the indicial equation of the recurrence `):
print(`Try:`):
print(`ShoreshABCa(N,1,1,1,1);`):

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



#########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:
option remember:
KAMA:=30:
 
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:=proc(n) local i:
lcm(seq(i,i=1..n)):
end:



#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) local gu: 
if type(a,integer) then
 RETURN(FAIL):
fi:

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

evalf(gu,20):

end:


##start guessing [K,r]

#GuessK1(L,m): Inputs a sequence of integers L, 
# finds a constant K of the form K=ithprime(i)^j, with i and j<=m   such that L[n]/K^n is an integer for all n. Try:
#GuessK1([seq((i^3+1)*25^(i),i=1..30)],10);
GuessK1:=proc(L,m) local K,i,lu,i1,p,K1,j:

K:=1:

for i from 1 to m do
p:=ithprime(i):

for j from m by -1  while max(seq(denom(L[i1]/(p^j)^i1),i1=1..nops(L)))>3000  to 0 do od:
 K:=K*p^j :
od:

K:

end:


#GuessK(L,m,r): inputs a sequence  of rational numbers, L, a positive integer m, and a positive integer r
#finds a rational number K such that L[i]*K^i*dn(r*i) are all integers. Try:
#GuessK([seq((2/3)^i/dn(5*i),i=1..20)],30,5);
GuessK:=proc(L,m,r) local L1,L1t,L1b,i,K1,K2,K,vu:

L1:=[]:

for i from 1 to nops(L) do
 L1:=[op(L1), L[i]*dn(r*i)]:
od:



L1t:=numer(L1):


L1b:=denom(L1):

K1:=GuessK1(L1t,m):


	
K2:=GuessK1(L1b,m):



K:=K2/K1:


	vu:= max(seq(denom(L[i]*dn(r*i)*K^i),i=1..nops(L)) ):


if vu>3000 then
#print(K, `did not work out`):
RETURN(FAIL):
fi:

K:

end:


#Nu1(a,b,K,r): if you have a sequence of diophantine approximations A(n)+B(n)*c such that
#(i)|A(n)+B(n)*c|<=Cb^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers
#outputs the implied  upper bound for the irrationality measure. Try:
#Nu1(op(ShoreshPi(N)[2]),25/16,10);
Nu1:=proc(a,b,K,r) local lu:
lu:=-(log(b)+r+log(K))/(log(a)+r+log(K)):

evalf(1+1/lu,20):
end:




#OpePiAB(A,B,n,N): the linear difference operator in backward notation for computing JF(4,2,5,6,10,A,B,2*n,x)
OpePiAB:=proc(A,B,n,N) local x: 
if A<=10 and B<=10 then
 DB(A,B,n,N):
elif A=8 and B=13 then
 Ope813(n,N):
else
 AZd1(JF(4,2,5,6,10,A,B,2*n,x),x,n,N): 
fi:

end:

#MakeDB(M,n,N): inputs a positive integer M, and outputs a database of all the operators for computing
#JF(4,2,5,6,10,A,B,2*n,x) for 1<=A,B<=N and gcd(A,B)=1. The output is a table such that T[A,B] is the desired operator.
#Try:
#MakeDB(2,n,N)
MakeDB:=proc(M,n,N) local T,A,B:

for A from 1 to M do
 for B from 1 to M do
  if igcd(A,B)=1 then
   T[A,B]:=OpePiAB(A,B,n,N):
  fi:
 od:
od:

op(T):

end:



#SeqFromRec(ope,n,N,INI,M): given a linear recurrence operator ope in n and the negative shift operator N^(-1) and
#a list of initial values (starting at n=0), outputs the first  M+1 terms of the sequence starting at n=0.
#Try:
#SeqFromRec(1/N+1/N^2,n,N,[0,1],10);
SeqFromRec:=proc(ope,n,N,INI,M) local T,n0,d,i:
d:=-ldegree(ope,N):

if nops(INI)<>d then
 RETURN(FAIL):
fi:

for n0 from 0 to d-1 do
 T[n0]:=INI[n0+1]:
od:

for  n0 from d to M do
 T[n0]:=add(subs(n=n0,coeff(ope,N,-i))*T[n0-i],i=1..d):
od:
[seq(T[n0],n0=0..M)]:
end:



#EmpDelABCa(A,B,C,a,M): the smallest empirical delta from M-10 to M for the Diophatine approximations implied by
#LABCa(A,B,C,a,n0). Try:
#EmpDelABCa(1,1,1,1,40);
EmpDelABCa:=proc(A,B,C,a,M) local gu,lu,i,lu1:

gu:=SeqPair(A,B,C,a,M):
if gu=FAIL then
 RETURN(FAIL):
fi:


lu:=min(seq(deltE(-gu[1][i]/gu[2][i],log((a+1)/a)), i=M-10..M)):

evalf(lu,20):
end:



#ShoreshABC(N,A,B,C,a): the indicial polynomial of the recurrence for Rukhadze's integral with A,B,C,a
#absolute values of the largest and smalleat roots of the indicial equation of the recurrence for the Salikhov Pi sequence
#linear recurrence operator with constant coefficients for Rukhdze's approximation
#Try:
#ShoreshABCa(N,1,1,1,1);
ShoreshABCa:=proc(N,A,B,C,a) local x,ope,lu,i,n:
option remember:

 ope:=AZd((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x,n,N)[1]:
ope:=lcoeff(ope,n):
ope:=normal(ope/igcd(seq(coeff(ope,N,i),i=0..degree(ope,N)))):

lu:=[solve(ope,N)]:
lu:=[seq(abs(lu[i]),i=1..nops(lu))]:
ope, evalf([max(abs(lu[1]),abs(lu[2]) ), min(abs(lu[1]),abs(lu[2]))],20):
end:


#Nu1(a,b,K,r): if you have a sequence of diophantine approximations A(n)+B(n)*c such that
#(i)|A(n)+B(n)*c|<=Cb^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers
#outputs the implied upper bound for the irrationality measure. Try:
#Nu1(op(ShoreshPi(N)[2]),25/16,10);
Nu1:=proc(a,b,K,r) local lu:
lu:=-(log(b)+r+log(K))/(log(a)+r+log(K)):

evalf(1+1/lu,20):
end:





#GuessKr(L,m): inputs a sequence  of rational numbers, L, a positive integer m,
#finds a positive integer r, and a  rational number K such that L[i]*K^i*dn(r*i) are all integers. Try:
#GuessKr([seq((2/3)^i/dn(5*i),i=1..20)],30);
GuessKr:=proc(L,m) local r,i,K:

if nops(L)<20 then
 print(`List too short, make it at least of length 20`):
 RETURN(FAIL):
fi:



if type(L[nops(L)],integer) then
 r:=0:
else

r:=round(ifactors(denom(L[nops(L)]))[2][-1][1]/nops(L));

for i from 15 to nops(L) do
  if denom(L[i])<>1 then

   if not member(round(ifactors(denom(L[i]))[2][-1][1]/i),{r,r-1}) then
 #   print(r, `did not work out`):
     RETURN(FAIL):
   fi:
 fi:
od:
fi:



K:=GuessK(L,m,r):




if K=FAIL then
# print(`r is`, r, `but could not find K`):
 RETURN(FAIL):
fi:

[r,K]:

end:


#GuessKrPair(P,m): Given a pair,P, of sequences of rational numbers and a positive integer m
#finds a non-negative integer r and a rational number K such that
#BOTH P[1][n]*dn(r*n)*K^(i*n) and P[1][n]*dn(r*n)*K^(i*n) are sequences of INTEGERS. 
#It tries out all the primes up to m and all the powers up to m.
#Type:

GuessKrPair:=proc(P,m) local gu1,gu2:

gu1:=GuessKr(P[1],m):

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

gu2:=GuessKr(P[2],m):

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

[max(gu1[1],gu2[1]), lcm(numer(gu1[2]), numer(gu2[2]))/igcd(denom(gu1[2]), denom(gu2[2]))]:

end:

#NorPair(P,r,K): Given a pair of sequences of rational numbers P, multiplies the n-th term by
#by lcm(1..r*n)*K^n. It also returns the list of gcds, let's call it G, and the list of log(G[i])/i and
#and the smallest and largest log(G[i])/i in the second half, in floating point. Try:
#NorPair(SeqPair(6,8,7,1,40),7,1);
NorPair:=proc(P,r,K) local P1,P2,vu,gu,i1:

if not (type(P,list) and nops(P)=2 and type(P[1],list) and type(P[2],list) and nops(P[1])=nops(P[2])) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if not (type(r,integer) and r>=0 and (type(K,integer) or type(K,fraction)) ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


P1:=[seq(P[1][i1]*dn(r*i1)*K^i1,i1=1..nops(P[1]))]:
P2:=[seq(P[2][i1]*dn(r*i1)*K^i1,i1=1..nops(P[2]))]:
vu:=max(op(denom(P1)),op(denom(P2))):

if vu>3000 then
 RETURN(FAIL):
fi:

P1:=vu*P1:
P2:=vu*P2:

gu:=[seq(igcd(P1[i1],P2[i1]),i1=1..nops(P1))]:

vu:=[seq(evalf(log(gu[i1])/i1, 20),i1=trunc(nops(P1)/2)..nops(P1))]:


[[P1,P2],gu,vu,min(op(vu)),max(op(vu))]:
end:






#Nu1Extra(a,b,K,r,K1): if you have a sequence of diophantine approximations A(n)+B(n)*c such that
#(i)|A(n)+B(n)*c|<=CONST*b^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers
#AND it looks like, empirically or rigorously that gcd(K^n*dn(r*n)*A(n), K^n*dn(r*n)*B(n)) are O(e^(K1*n))
#outputs the implied upper bound for the irrationality measure. Try:
#Nu1Extra(op(ShoreshPi(N)[2]),25/16,10,0);
Nu1Extra:=proc(a,b,K,r,K1) local lu:
lu:=-(log(b)+r+log(K)-K1)/(log(a)+r+log(K)-K1):

evalf(1+1/lu,20):
end:


#LABCa(A,B,C,a,n1): int((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x=0..1). Try:
#[seq(LABCa(1,1,1,1,n1),n1=0..10)];
#[seq(LABCa(6,8,7,1,n1),n1=0..10)];
LABCa:=proc(A,B,C,a,n1) local x:
int((x^A*(1-x)^B/(a+x)^C)^n1/(x+a),x=0..1):
end:

OpeABCa:=proc(A,B,C,a,n,N) local x:
option remember:
AZd((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x,n,N)[1]:
end:


OpeABCa1:=proc(A,B,C,a,n,N) local x:
option remember:
AZd1((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x,n,N):
end:

#SeqABCa(A,B,C,a,M): the first M+1 terms, starting at n=0, of the sequence LABCa(A,B,C,a,n1) via the Almkvist-Zeilberger recurrence. Try:
#SeqABCa(6,8,7,1,100);
SeqABCa:=proc(A,B,C,a,M) local ope,INI,n,N:
ope:=OpeABCa1(A,B,C,a,n,N):
INI:=[LABCa(A,B,C,a,0),LABCa(A,B,C,a,1)]:
SeqFromRec(ope,n,N,INI,M):
end:



#SeqPairOLD(A,B,C,a,M): outputs a pair of integer sequences of length M.
#The first one is the first M terms, starting at n=1, of the coefficient of 1 in
#LABCa(A,B,C,a,n1) and the second the coefficient of log(a+1). Try:
#SeqPairOLD(6,8,7,1,100);
SeqPairOLD:=proc(A,B,C,a,M) local gu,i,lu1,lu2,hal:


gu:=SeqABCa(A,B,C,a,M):

lu1:=[]:
lu2:=[]:
for i from 2 to M+1 do
 hal:=BU(gu[i],a):
 if hal=FAIL then
    RETURN(FAIL):
 else
  lu1:=[op(lu1),hal[1]]:
  lu2:=[op(lu2),hal[2]]:
 fi:
od: 
[lu1,lu2]:

end:




#SeqPair(A,B,C,a,M): outputs a pair of integer sequences of length M.
#The first one is the first M terms, starting at n=1, of the coefficient of 1 in
#LABCa(A,B,C,a,n1) and the second the coefficient of log(a+1). Try:
#SeqPair(6,8,7,1,100);
SeqPair:=proc(A,B,C,a,M) local X,ope,INI,n,N,i,lu0,gu,gu0:

gu0:=SeqABCa(A,B,C,a,1)[2]:

lu0:=BU(gu0,a):

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

ope:=OpeABCa1(A,B,C,a,n,N):
INI:=[X,lu0[1]+X*lu0[2]]:
gu:=SeqFromRec(ope,n,N,INI,M):

[[seq(coeff(gu[i],X,0),i=2..M+1)], [seq(coeff(gu[i],X,1),i=2..M+1)]]:

end:


#BestABCa(a,K,M): Explores the empirical delta for the terms between n=M-10 and n=M, in the rational approximations in the Rukhazde integral
#LABCa(A,B,C,a,n1) for 1<=A,B,C<=K whose delta is >cu, Try
#BestABCa(1,3,500);
BestABCa:=proc(a,K,M) local cu,gu,A1,B1,C1,A,B,C,x,n,aluf, si,halev,t0:
t0:=time():

gu:=(x^A*(1-x)^B/(a+x)^C)^n/(x+a):

print(`The best choice of the Rukhazde integral `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In  Rukhazde's beautiful article `):
print(``):
print(subs({A=6,B=8,C=7},gu)):
print(``):
print(`This inspired us to search for variations where the exponents 6,8,7 are replaced by other choices`):
print(`A,B,C relatively prime between 1 and`, M):
print(``):
print(gu):
print(``):

print(`Since a rigorous investigation (even by computer) for each case is painful, it makes sense to first find the empirical deltas in `):
print(``):
print(`abs(Pi-an/bn)<=CONSTANT/bn^(1+delta) `):
print(``):
print(` for fairly large n. We call this delta the empirical delta, and take, in each case, the smallest for n between`, M-10, `and ` , M):
print(``):
print(`This will give us an indication of what values of A and B are promising and worth pursuing. `):
print(``):
print(`Thanks to the recurrences produced, in each case, by the amazing Almkvist-Zeilberger algorithm, it is very easy to crank-out many terms`):
print(``):
print(`of these integrals, much faster then a direct computation. `):
print(``):

print(`For the sake of completeness we will list ALL the empirical deltas that exceed the Alladi-Robinson choice of A=1,B=1,C=1`):
print(``):


aluf:=[1,1,1]:

si:=EmpDelABCa(1,1,1,1,M):

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

cu:=si:

for A1 from 1 to K do
 for B1 from 1 to K do
  for C1 from 1 to K do
  if igcd(A1,B1,C1)=1 then

  halev:=EmpDelABCa(A1,B1,C1,a,M):
     if halev<>FAIL and halev>cu then
  print(`(A,B,C)=`, (A1,B1,C1), `EmpDel=`, halev):
  if halev>si then
    si:=halev:
     aluf:=[A1,B1,C1]:
  fi:
 fi:
 fi:
 od:
od:
od:


print(`The highest Empirical delta is `, si, `achieved by the case A=`, aluf[1], `and B= `, aluf[2], `and C=`, aluf[3] ):
print(``):
print(`this leads to an irrationality measure around`, evalf(1+1/si,20)):


print(``):
print(`Hence it is worthwhile to try to investigate this case of`, aluf):
print(``):
print(``):
print(`-----------------------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to generate. `):
print(``):

end:


#BU(f,a): If f is A+B*log((a+1)/a)) with A and B rational functions returns [A,B]
#Try
#BU(11/19*log(5/4)-13/109);
BU:=proc(f,a) local f1,i,A,B:
f1:=expand(f):

if type(f1,function) then
 if f1=log((a+1)/a) then
   RETURN([0,1]):
 fi:
fi:


if type(f1,`*`) then
if type(op(2,f1),function) then
 if op(2,f1)=log((a+1)/a) then
   RETURN([0,op(1,f1)]):
 fi:
fi:
fi:

if not type(f1,`+`) then
 RETURN(FAIL):
fi:

for i from 1 to nops(f1) do
 if (type(op(i,f1),fraction) or  type(op(i,f1),integer) ) then
  A:=op(i,f1):
  B:=(f1-A)/log((a+1)/a):
  B:=simplify(B):
  if not (type(B,integer) or  type(B,fraction)) then
   RETURN(FAIL):
  fi:

  if simplify(f-A-B*(log(a+1)-log(a)))<>0 then
     RETURN(FAIL):
   else
     RETURN([A,B]):
    fi:
 fi:
od:

   A:=0:
  B:=simplify((f1-A)/log((a+1)/a)):
  if simplify(f-A-B*(log(a+1)-log(a)))<>0 then
     RETURN(FAIL):
  fi:
 [A,B]:
end:


#InfoABCa(A,B,C,a,x,n,N,M): 
#InfoABCa(6,8,7,1,n,N,400);
InfoABCa:=proc(A,B,C,a,x,n,N,M) local ope,ope1,gu,lu0,lu1,lu2,lu3,lu4,lu5,kak:
option remember:

if not (type(A,integer) and type(B,integer) and type(C,integer) and igcd(A,B,C)=1 and type(n,symbol) and type(N,symbol) and type(M,integer)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=[Int((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x=0..1)]:
lu0:=EmpDelABCa(A,B,C,a,M):
if lu0=FAIL then
 RETURN(FAIL):
fi:

lu0:=evalf(1+1/lu0,20):
gu:=[op(gu),lu0]:


ope:=AZd((x^A*(1-x)^B/(a+x)^C)^n/(x+a),x,n,N)[1]:

ope1:=lcoeff(ope,n):
ope1:=numer(ope1/coeff(ope1,N,0)):
gu:=[op(gu),ope,ope1]:


lu1:=ShoreshABCa(N,A,B,C,a)[2]:
gu:=[op(gu),lu1]:

 
kak:=SeqPair(A,B,C,a,40):

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

lu2:=GuessKrPair(kak,20);

if lu2=FAIL then
 RETURN(gu):
else
gu:=[op(gu),lu2]:
fi:



lu3:=Nu1Extra(lu1[1],lu1[2],lu2[2],lu2[1],0):

gu:=[op(gu),lu3]:



lu4:=NorPair(SeqPair(A,B,C,a,M),lu2[1],lu2[2])[4];

gu:=[op(gu),lu4]:


lu5:=Nu1Extra(lu1[1],lu1[2],lu2[2],lu2[1],lu4):

gu:=[op(gu),lu5]:

gu:

end:

#PaperABCa(A,B,C,a,M):
PaperABCa:=proc(A,B,C,a,M) local x,n,N,gu,ope,ope1,F,i1,alef,bet,alef1,bet1,D,C1,D1,r,K,lu,K1,t0,n0:
t0:=time():

gu:=InfoABCa(A,B,C,a,x,n,N,M):

if gu=FAIL or nops(gu)<9 then
 print(`Sorry, no paper`):
 RETURN(FAIL):
fi:

if gu[9]>0 then

print(`A sequence of Diophantine Approximations to`, log(a+1)-log(a), ` that implies its irrationality`):
print(`With a fully rigorous crude upper bound on the irrationality measure of`, gu[7] ):
print(`and a not-yet-rigorous smaller upper bound around`, gu[9]):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
else
print(`A sequence of Diophantine Approximations to Pi that Disappointingly does not imply its irrationality`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
fi:

print(`Consider the following Rukhazde-type integral, let's call it E(n) `):
print(``):
print(gu[1]):
print(``):
print(`It is readily seen that E(n) can be written as`):
print(``):
print(`E(n)=C(n) + D(n)*log((a+1)/a) `):
print(``):
print(`For certain sequences of RATIONAL numbers C(n), D(n)`):
print(``):
ope:=gu[3]:
print(`It follows from the amazing Almkvist-Zeilberger algorithm that E(n), and hence C(n) and D(n) satisfy the following`):
print(` second-order linear recurrence equation with polynomial coefficients. `):
print(``):
print(add(coeff(ope,N,i1)*F(n+i1),i1=0..degree(ope,N))=0):
print(``):
print(`and in Maple format `):
print(``):
lprint(add(coeff(ope,N,i1)*F(n+i1),i1=0..degree(ope,N))=0):
print(``):
print(`The initial conditions are as follows`):
print(``):

print(seq( E(n0)=LABCa(A,B,C,a,n0),n0=0..1)):
print(``):
print(`and in Maple notation`):
print(``):
print(seq( E(n0)=LABCa(A,B,C,a,n0),n0=0..1)):
print(``):

print(`Using this recurrence, we can compute many values, and use them to ESTIMATE the implied bound irrationality measure.`):
print(`Using the values from  n=`, M-10, ` to  M=`, M, `and taking the largest value yields the following estimate`):
print(``):
lprint(gu[2]):
print(``):
print(`We can use the Poincare lemma to find the asymptotic behavior of the exponentially growing C(n), D(n) and`):
print(`the exponentially decaying E(n).`):
print(``):
ope1:=gu[4]:


print(`The constant-coefficient approximation to the above recurrence has indicial polynomial `):
print(``):
print(ope1):
print(``):
print(`and in Maple format `):
print(``):
lprint(ope1):
print(``):
print(`whose largest root let's call it alef, is`, gu[5][1], `and absolute value of the  smaller root, let's call it bet is`, gu[5][2]):
print(``):
print(`It follows that, up to polynomial correction that will get out in the wash,  C(n) and D(n) are OMEGA(alef^n) and E(n) is O(bet^n)`):
print(``):
print(`We now need a divisibility lemma that is left to the reader`):
print(``):
print(`Let d(n)=lcm(1..n) `):
r:=gu[6][1]:
K:=gu[6][2]:


print(`Lemma: `):
print(C1(n)=d(r*n)*K^n*C(n),  D1(n)=d(r*n)*K^n*D(n)):
print(`are integers`):
print(`Defining `):
print(``):
print(E1(n)=d(r*n)*K^n*E(n)):
print(``):
print(`We have `):
print(``):
print(E1(n)=C1(n)+ log((a+1)/a)*D1(n)):
print(``):
print(`where C1(n) and D1(n) are INTEGERS `):
print(``):


print(`Since, famously, d(n)=OMEGA(e^n), we have `):
print(``):
print(`C1(n) and D1(n) are  OMEGA of`):
print(exp(r*n)*K^n*alef^n):
print(``):
print(`and E1(n) is O of `):
print(``):
print(exp(r*n)*K^n*bet^n):
print(``):
alef1:=gu[5][1]:
bet1:=gu[5][2]:


print(`It follows that  E1(n)=max(C1(n),D1(n))^(-delta) where delta equals `):
print(``):
lu:=-(log(bet)+r+log(K))/(log(alef)+r+log(K)):

print(lu):
print(`recall that alef is `, alef1, `and bet is `, bet1):
lu:=evalf(subs({bet=bet1,alef=alef1},lu),20):
print(`and delta is `):
print(``):
lprint(lu):
print(``):
print(`implying the rigorous, but crude upper bound for the irrationality measure, 1+1/delta that equals `):
print(``):
lprint(evalf(1+1/lu,20)):
print(``):
print(`confirming the first part of the statement in the title `):
print(``):
print(`We now need a harder lemma, that we confess that we can't do, but you the reader may be able to `):
print(``):
print(`Harder Lemma: There exists a rigorously proved constant K1 such that `):
print(``):
print(`gcd(C1(n),D1(n))=O(exp(K1*n)) `):
print(``):
print(`Once we find such a constant K1 the delta can be improved to `):
print(``):
lu:=-(log(bet)+r+log(K)-K1)/(log(alef)+r+log(K)-K1):
print(``):
print(`By looking at the smallest value of log(gcd(C1(n),D1(n))/n for n between n=`, M/2, `and `, M, `we estimate that K1 can be taken to be`):
print(``):
lprint(gu[8]):
print(``):
print(`and the better delta is `):
print(``):
lu:=subs({alef=alef1,bet=bet1,K1=gu[8]},lu):
lu:=evalf(lu,20):
print(``):
lprint(lu):
print(``):
print(`that implies an irrationality measure, 1+1/delta that equals `):
print(``):
lprint(evalf(1+1/lu,20)):
print(``):
print(`-----------------------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to generate. `):
print(``):

end:





#BU(f,a): If f is A+B*log((a+1)/a)) with A and B rational functions returns [A,B]
#Try
#BU(11/19*log(5/4)-13/109);
BU:=proc(f,a) local f1,i,A,B:
f1:=expand(f):

if type(f1,function) then
 if f1=log((a+1)/a) then
   RETURN([0,1]):
 fi:
fi:


if type(f1,`*`) then
if type(op(2,f1),function) then
 if op(2,f1)=log((a+1)/a) then
   RETURN([0,op(1,f1)]):
 fi:
fi:
fi:

if not type(f1,`+`) then
 RETURN(FAIL):
fi:

for i from 1 to nops(f1) do
 if (type(op(i,f1),fraction) or  type(op(i,f1),integer) ) then
  A:=op(i,f1):
  B:=(f1-A)/log((a+1)/a):
  B:=simplify(B):
  if not (type(B,integer) or  type(B,fraction)) then
   RETURN(FAIL):
  fi:

  if simplify(f-A-B*(log(a+1)-log(a)))<>0 then
     RETURN(FAIL):
   else
     RETURN([A,B]):
    fi:
 fi:
od:

   A:=0:
  B:=simplify((f1-A)/log((a+1)/a)):
  if simplify(f-A-B*(log(a+1)-log(a)))<>0 then
     RETURN(FAIL):
  fi:
 [A,B]:
end:











#BestABCaG(a,K,M): Explores the empirical delta for the terms between n=M-10 and n=M, in the rational approximations in the Rukhazde integral
#LABCaG(A,B,C,a,n1) for 1<=A,B,C<=K where A and B are between C-G and C+G whose delta is  larger than the one for (1,1,1), Try
#BestABCaG(1,3,500,1);
BestABCaG:=proc(a,K,M,G) local cu,gu,A1,B1,C1,A,B,C,x,n,aluf, si,halev,t0:
t0:=time():

gu:=(x^A*(1-x)^B/(a+x)^C)^n/(x+a):

print(`The best choice of the Rukhazde integral `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`consider the integral `):
print(gu):
print(``):
print(`When A=B=C=1 it is the Alladi-Robinson integral.  Can we do better?`):
print(``):

print(`Since a rigorous investigation (even by computer) for each case is painful, it makes sense to first find the empirical deltas in `):
print(``):
print(`abs(Pi-an/bn)<=CONSTANT/bn^(1+delta) `):
print(``):
print(` for fairly large n. We call this delta the empirical delta, and take, in each case, the smallest for n between`, M-10, `and ` , M):
print(``):
print(`This will give us an indication of what values of A and B are promising and worth pursuing. `):
print(``):
print(`Thanks to the recurrences produced, in each case, by the amazing Almkvist-Zeilberger algorithm, it is very easy to crank-out many terms`):
print(``):
print(`of these integrals, much faster then a direct computation. `):
print(``):

print(`For the sake of completeness we will list ALL the empirical deltas that exceed the Alladi-Robinson choice of A=1,B=1,C=1`):
print(``):


aluf:=[1,1,1]:

si:=EmpDelABCa(1,1,1,a,M):

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

print(`Our base line is the value for A=1,B=1,C=1 that is`, si):

cu:=si:

  for C1 from 1 to K do

   for A1 from max(1,C1-G) to C1+G  do
     for B1 from  max(1,C1-G) to C1+G  do

  if igcd(A1,B1,C1)=1 then

  halev:=EmpDelABCa(A1,B1,C1,a,M):
     if halev<>FAIL and halev>cu then
  print(`(A,B,C)=`, (A1,B1,C1), `EmpDel=`, halev):
  if halev>si then
    si:=halev:
     aluf:=[A1,B1,C1]:
  fi:
 fi:
 fi:
 od:
od:
od:


print(`The highest Empirical delta is `, si, `achieved by the case A=`, aluf[1], `and B= `, aluf[2], `and C=`, aluf[3] ):
print(``):
print(`this leads to an irrationality measure around`, evalf(1+1/si,20)):


print(``):
print(`Hence it is worthwhile to try to investigate this case of`, aluf):
print(``):
print(``):
print(`-----------------------------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `to generate. `):
print(``):

end: