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

Digits:=200:

print(` Written: Dec. 2020- Jan. 2021 `):
print():
print(`This is GenBeukersLog.txtG, A Maple package for efficient computation of generalized Beukers integrals`):
print(`that generalize the Alladi-Robinson log integrals `):

print(` It experiments with the Generalized Alladi-Robinson integral`):
print(`Int(x^(a)*(1-x)^(b)*(x*(1-x)^n)/(1+c*x)^(n+d+1),x=0..1):`):
print(``):
print(`where the case c=1 d=0, correspond to the former case`):

print(` It is one of the Maple packages that accompany   the article `):
print(`Tweaking the Beukers Integrals In Search of More Miraculous Irrationality Proofs A La Apery`):
print(``):
print(`avaliable from the authors' websites and from arxiv.org `):
print(``):

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



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


ezraGen:=proc()
if args=NULL then
 print(`The Generarilzed procedures are:  AnBnG, CnDnG, FindRelG, GLCg,  Hopefuls2, Hopefuls2C, Hopefuls2v, IntGLg, OpeLg `):
else
 ezra(args):
fi:

end:

ezraPrimes:=proc()
if args=NULL then
 print(`The Primes procedures are: AsyPpG, EvalPpG, Fs, Lc,PpG , PrimesF, PrintPpG, Seg1 `):
else
 ezra(args):
fi:

end:

ezra1:=proc()
if args=NULL then
 print(`The supsporting procedures are: Hopefuls1,Hopefuls1v, FindRel, IntGLsummand, IntGLsummand, IntGL2F1, Lcons, IR, MyIDs, Search, SeqFromRec`):
else
 ezra(args):
fi:

end:

ezra:=proc()
if args=NULL then

print(`The MAIN procedures are:`) :
print(` AnBn, AppxSeq, CnDn, delt, deltSeq, GLC, Info1, IntGL, OpeL, TheoremL  `):


elif nargs=1 and args[1]=AnBn then
print(`AnBn(a,b,c,K): The first K terms in the sequences An and Bn in the expression`):
print(`IntGL(a,b,c,n)=B(n)*C-A(n)`):
print(`where C=IntGL(a,b,c,0); `):
print(`It also returns the first K terms in the floating point of the sequence of errors. Try: `):
print(`AnBn(0,0,1,50);`):


elif nargs=1 and args[1]=AnBnG then
print(`AnBnG(a,b,c,d,K): The first K terms in the sequences An and Bn in the expression`):
print(`IntGLg(a,b,c,d,n)=B(n)*C-A(n)`):
print(`where C=IntGLg(a,b,c,d,0); `):
print(`It also returns the first K terms in the floating point of the sequence of errors. Try: `):
print(`AnBnG(0,0,1,0,50);`):


elif nargs=1 and args[1]=AppxSeq then
print(`AppxSeq(a,b,c,K): The approximating sequence for  IntL(a,b,c,0); `):
print(`Try: `):
print(`AppxSeq(0,0,1,40);`):

elif nargs=1 and args[1]=AsyPpG then
print(`AsyPpG(PA,n1): Inputs a GpP expression in the form [Cn,Dn] where Cn a, Dn are the partial PpG expressions for the sequences Cn and Dn (see paper)`):
print(`outputs nu sich that the INTEGERating sequence Dn/Cn is asympotic to exp(nu*n). `):
print(`For GLC(0,0,1) (alias log(2))`):
print(`AsyPpG([ [[1,1],[1,0],[]],    [[1, 1], [1, 1],[]  ] ] );`):
print(``):
print(`For GLC(1/3,1/3,1/5) `):
print(`AsyPpG([ [[5,1],[1,0],[  [1/3,2/3,0,3,{1},3   ]  ]],    [[3, 3/2], [1, 0],[[0,2/3,0,3,{2},3   ]  ] ]] );`):
print(``):
print(`For GLC(1/2,1/2,1) `):
print(`AsyPpG([ [[1,0],[1,0],[  [1/2,1,0,2,{0},1   ]  ]],    [[2, 2], [1, 0],[  ] ]] );`):

print(``):
print(`For GLC(1/3,1/3,1) `):
print(`AsyPpG([ [[1,0],[1,0],[  [1/3,2/3,0,3,{1},3   ]  ]],    [[3, 3/4], [1, 0],[[0,2/3,0,3,{2},3]   ] ]] );`):


print(``):
print(`For GLC(1/4,1/4,1) `):
print(`AsyPpG([ [[1,0],[1,0],[  [1/4,1/2,0,4,{1},4   ]  ]],    [[2, 3], [1, 0],[[0,3/4,0,2,{3},4]   ] ]] );`):

elif nargs=1 and args[1]=CnDn then
print(`CnDn(a,b,c,K): The first K terms in the sequences gcd(numer(An[i]),number(Bn[i]))`):
print(`lcm(denom(An[i]),denom(Bn[i]))`):
print(`It also returns the last 20 terms of the logs of the normailizing sequence`):
print(`as well as the last 20 terms of the implied deltas`):
print(`Try:`):
print(`CnDn(0,0,1,100);`):


elif nargs=1 and args[1]=CnDnG then
print(`CnDnG(a,b,c,d,K): The first K terms in the sequences gcd(numer(An[i]),number(Bn[i]))`):
print(`lcm(denom(An[i]),denom(Bn[i]))`):
print(`It also returns the last 20 terms of the logs of the normailizing sequence`):
print(`as well as the last 20 terms of the implied deltas`):
print(`Try:`):
print(`CnDnG(0,0,1,0,100);`):

elif nargs=1 and args[1]=EvalPpG then
print(`EvalPpG(PA,n1): evaluates a   PpG [Cn,Dn] where Cn and Dn are partial PpG objects`):
print(`For GLC(0,0,1) (alias log(2))`):
print(`EvalPpG( [ [[1,1],[1,0],[]],    [[1, 1], [1, 1],[] ] ] ,1000  );`):

elif nargs=1 and args[1]=FindRel then
print(`FindRel(a,b,c): A  relationship`):
print(`between IntGL(a,b,c,0) and  IntGL(a,b,c,1) .  The output is of the form [[a0,a1],c] , that means that`):
print(`a0*IntGL(a,b,c,0)+a1*IntGL(a,b,c,1)=c. Try:`):
print(`FindRel(0,0,1);`):


elif nargs=1 and args[1]=FindRelG then
print(`FindRelG(a,b,c,d): A  relationship`):
print(`between IntGLg(a,b,c,d,0) and  IntGLg(a,b,c,d,1) .  The output is of the form [[a0,a1],c] , that means that`):
print(`a0*IntGL(a,b,c,0)+a1*IntGL(a,b,c,1)=c. Try:`):
print(`FindRelG(0,0,0,1);`):


elif nargs=1 and args[1]=Fs then
print(`Fs(S,n): Given a set S and a positive integer n finds sort([seq(evalf(frac(n/p),10),p in S)]);`):


elif nargs=1 and args[1]=GLC then
print(`GLC(a,b,c): A floating-point approximation to the constant`):
print(`IntGL(a,b,c,0) supposed to be good for the number of digits. Try:`):
print(`GLC(0,0,1);`):


elif nargs=1 and args[1]=GLCg then
print(`GLCg(a,b,c,d): A floating-point approximation to the constant`):
print(`IntGLg(a,b,c,d,0) supposed to be good for the number of digits. Try:`):
print(`GLCd(0,0,1,0);`):



elif nargs=1 and args[1]=Hopefuls1 then
print(`Hopefuls1(): A list of length 9 where the entries are lists of hopeful triples [a,b,c], with denominators 1,2,3,4,5,6,7,8,9, respectively,`):
print(`and a,b,c  are positive`):
print(`It is those triples where CnDn(a,b,c,100)[3] are all positive. In other words, it is the pre-computed output of`):
print(`Search(N,100) for N=1..9 `):
print(`Try: `):
print(` Hopefuls1(); `):

elif nargs=1 and args[1]=Hopefuls1v then
print(`Hopefuls1v(K): Lots of constants that probably have an irrationality proof, together for each the relevant information. Try:`):
print(`Hopefuls1v(2000):`):


elif nargs=1 and args[1]=Hopefuls2 then
print(`Hopefuls2(): A list of length 5 where the entries are lists of hopeful of 4-tuples [a,b,c,d], with denominators 1,2,3,4,5 respectively,`):
print(`and a,b,c,d  are positive`):
print(`It is those triples where CnDnG(a,b,c,d,1000)[3] are all positive. In other words, it is the pre-computed output of`):
print(`Try: `):
print(` Hopefuls2(); `):


elif nargs=1 and args[1]=Hopefuls2C then
print(`HopefulsC(r): A list of length 5 where the entries are lists of hopeful of 4-tuples [a,b,c,d], with denominators 1,2,3,4,5 respectively,`):
print(`and a,b,c,d  are positive`):
print(`It is those triples where CnDnG(a,b,c,d,1000)[3] are all positive. In other words, it is the pre-computed output of`):
print(`Try: `):
print(` Hopefuls2C(r); `):

elif nargs=1 and args[1]=Hopefuls2v then
print(`Hopefuls2v: A description of generalized hopefuls. Try:`):
print(`Hopefuls2v():`):

elif nargs=1 and args[1]=Info1 then
print(`Info1(a,b,c,K): inputs a triple (a,b,c) and outputs a list consisting of`):
print(`(i) the constant as a 2F1, (ii) the floating point to 100 digits,`):
print(`(iii) the recurrence operator annihilating An and Bn (iv) the pair of initial conditions for An and Bn`):
print(`i.e. [[A(0),A(1)],[B(0),B(1)]] (v) the estimated value of the limit of the log of the INTEGERATING functor divided by n (vi) the`):
print(`estimated implied irrationality measure using K values. Try:`):
print(`Info1(0,0,1,2000);`):

elif nargs=1 and args[1]=IntGL then
print(`IntGL(a,b,c,n1):`):
print(`int( x^(a)*(1-x)^(b)*(x*(1-x))^n1/(1+c*x)^(n1+1); Try:`) :
print(`IntGL(0,0,1,0);`):


elif nargs=1 and args[1]=IntGLg then
print(`IntGLg(a,b,c,d,n1):`):
print(`int( x^(a)*(1-x)^(b)*(x*(1-x))^n1/(1+c*x)^(n1+1+d);  divided by int(x^a*(1-x)^b/(1+c*x)^d,x=0..1). Try:`) :
print(`IntGLg(0,0,1,0,0);`):


elif nargs=1 and args[1]=IntGLsummand then
print(`IntGLsummand(a,b,c,n): The summand in the series representation for IntGL(a,b,c) in terms of factrials. Try:`):
print(`IntGLsummand(0,0,1,n);`):

elif nargs=1 and args[1]=IntGLsummandRF then
print(`IntGLsummandRF(a,b,c,n,RF): The summand in the series representation for IntGL(a,b,c) in terms of the raising-factorial`):
print(`denoted by RF. Try`):
print(`IntGLsummandRF(0,0,1,n,RF);`):


elif nargs=1 and args[1]=IntGL2F1 then
print(`IntGL2F1(a,b,c): IntGL2F1(a,b,c): IntGL(a,b,c,0) expressed as a 2F1. Try:`):
print(`IntGL2F1(0,0,1);`):

elif nargs=1 and args[1]=IR then
print(`IR(L): IR(L): same as IntegerRelations[PSLQ](L). Try`):
print(` phi:=evalf((sqrt(5)+1)/2): IR([1,phi,phi^2]); `):

elif nargs=1 and args[1]=Lc then
print(`Lc(n): the lcm of 1, ..., n. Try: Lc(1000);`):


elif nargs=1 and args[1]=Lcons then
print(`Lcons(c): The pair of numbers [beta1,beta2] such that An,Bn=Omega(beta1)^b, I(n)=A(n)-c*B(n)= Omega(beta2^n).`):


elif nargs=1 and args[1]=MyIDs then
print(`MyIDs(C,F,F0,N): Given a constant C in decimals and another constant`):
print(`let's call it F, whose floating-point if F0`):
print(` and a positive integer N`):
print(`tries to express C as (a*F+b)/(c*F+d) for a,b,c,d from -N to N using PSLQ`):
print(`MyIDs(evalf((log(2)-2)/(2*log(2)+3)),log(2),evalf(log(2)),100);`):

elif nargs=1 and args[1]=OpeL then
print(`OpeL(a,b,c,n,N):`):
print(`The second-order linear recurrence operator annihinationg Int(x^a*(1-x)^b*(x*(1-x))^n/(1+c*x)^(n+1),x=0..1): Try`):
print(`OpeL(0,0,1,n,N);`):

elif nargs=1 and args[1]=OpeLg then
print(`OpeLg(a,b,c,d,n,N): The operator in n and the shift operator N, annihaling int(x^a*(1-x)^b/(1+c*x)^d (x*(1-x)/(1+c*d))^n,x=0..1):`):
print(`Try:`):
print(`OpeLg(0,0,1,0,n,N);`):

elif nargs=1 and args[1]=PpG then
print(`PpG(a,b,n,m1,m2,C,M): Given rational numbers a and b between 0 and 1 and positive integers m1,m2 and n`):
print(`and a subset C in {1,...M-1} and an integer M`):
print(`it is the product of p between m1 and m2  such that frac{n/p} is between a and b`):
print(`and p mod M is in C`):

elif nargs=1 and args[1]=PrimesF then
print(`PrimesF(N,n,k): inputs a HUGE integer N that is a product of small primes, and a much smaller integer n`):
print(`that generated N (via some process), outputs the list of length k, whose i-th entry is the list`):
print(`of primes>=sqrt(n) that show up with exponent i, followed by the list of primes that`):
print(`did not show up at all. Try:`):
print(`lu:=CnDn(1/3,1/3,1/15,1000)[1][2][-1]: PrimesF(lu,1000,6);`):


elif nargs=1 and args[1]=PrintPpG then
print(`PrintPpG(PA,n,L1,L2): P is the  same input as the P EvalPpG(P,n1) (q.v.) but n,L1 and L2 are symbolcs`):
print(`where n stands for n,  L1(n) stands for lcm(1...n) and L2(alpha,beta,a,b,C,M,n) stands for the`):
print(`product of all primes between a*n and b*n such that the fractional part of n/p is between alpha and beta`):
print(`and p mod M belongs to the set C where C is a subset of the set of nonneg. integers less than M relatively prime to it.`):
print(`If no modularity condition is present than M=1 and C={0}. Try:`):


elif nargs=1 and args[1]=Search then
print(`Search(N,K): All the triples [a,b,c]=[a'/N, b'/N, c'/N] with a', b' between 0 and N-1 and`):
print(`c' between 1 and N such that the deltas for CnDn(a,b,c) seems to be positive.`):
print(`Try:`):
print(`Search(2,300);`):

elif nargs=1 and args[1]=Seg1 then
print(`Seg1(L,a,b): Given a  list L and integers a<b, extracts the segment between a and b. Try:`):
print(`Seg1([seq(i,i=1..100)],40,60);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,ini,L): Given an (ordinary) recurrence operator`):
print(`ope(n,N) in the variable n and the shift operator N (Nf(n):=f(n+1))`):
print(`and given initial values [ini0,ini1,..,ini_{ORD-1}],computes`):
print(`the first L+1 terms of the sequence a(n) satisfying`):
print(`ope(n,N)a(n)=0 and a[i]=ini[i] for i=0,...,ORD-1`):
print(`For example, try: SeqFromRec(N^2-N-1,n,N,[0,1],10);`):

elif nargs=1 and args[1]=TheoremL then
print(`TheoremL(a,b,c,K,P): Inputs rational numbers a,b,c, and a large positive integer K and outputs`):
print(`a theorem regarding the constant`):
print(`print(Int(x^(a)*(1-x)^(b)/(1+c*x),x=0..1)/Int(x^(a)*(1-x)^(b),x=0..1).`):
print(`Either a suggested proof of irrationality or a way to compute it exponentially  fast.`):
print(`P is 0 if no INTEGERating factor is known, otherwise it is a proposed one.`):
print(`Try: `):
print(`For log(2) try: `):
print(`TheoremL(0,0,1,2000,  [ [[1,1],[1,0],[]],    [[1, 1], [1, 1],[]  ] ] ):`):
print(``):
print(`For GLC(1/3,1/3,1/5) `):
print(`TheoremL(1/3,1/3,1/5,2000,[ [[5,1],[1,0],[  [1/3,2/3,0,3,{1},3   ]  ]],    [[3, 3/2], [1, 0],[[0,2/3,0,3,{2},3   ]  ] ]] );`):
print(``):
print(`For GLC(1/2,1/2,1) `):
print(`TheoremL(1/2,1/2,1,2000,[ [[1,0],[1,0],[  [1/2,1,0,2,{0},1   ]  ]],    [[2, 2], [1, 0],[  ] ]] );`):

print(``):
print(`For GLC(1/3,1/3,1) `):
print(`TheoremL(1/3,1/3,1,2000,[ [[1,0],[1,0],[  [1/3,2/3,0,3,{1},3   ]  ]],    [[3, 3/4], [1, 0],[[0,2/3,0,3,{2},3]   ] ]] );`):


print(``):
print(`For GLC(1/4,1/4,1) `):
print(`TheoremL(1/4,1/4,1,2000,[ [[1,0],[1,0],[  [1/4,1/2,0,4,{1},4   ]  ]],    [[2, 3], [1, 0],[[0,3/4,0,2,{3},4]   ] ]] );`):

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

fi:


end:



#SeqFromRec(ope,n,N,ini,L): Given an (ordinary) recurrence operator
#ope(n,N) in the variable n and the shift operator N (Nf(n):=f(n+1))
#and given initial values [ini0,ini1,..,ini_{ORD-1}],computes
#the first L+1 terms of the sequence a(n) satisfying
#ope(n,N)a(n)=0 and a[i]=ini[i] for i=0,...,ORD-1
#For example, try: SeqFromRec(N^2-N-1,n,N,[0,1],10);

SeqFromRec:=proc(ope,n,N,ini,L)
local gu,lu,j,ORD,n0,i:

ORD:=degree(ope,N):
if ORD<>nops(ini) then
 ERROR(`The length of`,ini, `should be `, degree(ope,N)):
fi:
gu:=ini:
for i from ORD to L do
lu:=0:
n0:=i-ORD:
lu:=0:
for j from 0 to ORD-1 do
lu:=lu+subs(n=n0,coeff(ope,N,j))*gu[n0+j+1]:
od:
lu:=-lu/subs(n=n0,coeff(ope,N,ORD)):
gu:=[op(gu),lu]:
od:
gu:
end:

#IntGL(a,b,c,n1);
#intGL( x^(a)*(1-x)^(b)*(x*(1-x))^n/(1+c*x)^(n+1);
#IntGL(0,0,1,0);
IntGL:=proc(a,b,c,n) local x,y,lu,mu:
lu:=int(x^a*(1-x)^b*x^n*(1-x)^n/(1+c*x)^(n+1),x=0..1):
mu:=int(x^a*(1-x)^b,x=0..1):
lu/mu:
end:


#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): evalf(-log(abs(c-a))/log(denom(a)))-1: end:




#FindRel(a,b,c): A  relationship
#between IntGL(a,b,c,0) and  IntGL(a,b,c,1) .  The output is of the form [[a0,a1],c] , that means that
#a0*IntGL(a,b,c,0)+a1*IntGL(a,b,c,1)=c. Try:
#FindRel(0,0,1);
FindRel:=proc(a,b,c) local x,c0,c1,gu,eq,var,lu,A,B,C,gu1,gu2,hal,ka,i:

option remember:

lu:=x^(a)*(1-x)^(b):

gu1:=normal(diff((c0+c1*x)*x*(1-x)*lu/(1+c*x),x)):

gu1:=normal(gu1/lu):

gu2:=normal(A/(1+c*x) +B*x*(1-x)/(1+c*x)^2-C):

gu:=normal(gu1-gu2):

gu:=numer(gu):

var:={c0,c1,A,B,C}:

eq:={seq(coeff(gu,x,i),i=0..degree(gu,x))}:

var:=solve(eq,var):

hal:=subs(var,[A,B,C]):

if hal[1]=0 then
 RETURN(FAIL):
fi:

hal:=[1,normal(hal[2]/hal[1]),normal(hal[3]/hal[1])]:

hal:

ka:=lcm(seq(denom(hal[i]),i=1..nops(hal))):

hal:=ka*hal:

[[hal[1],hal[2]],hal[3]]:

end:


#OpeL(a,b,c,n,N):
#The second-order linear recurrence operator annihinationg Int(x^a*(1-x)^b*(x*(1-x))^n/(1+c*x)^(n+1),x=0..1): Try
#OpeL(0,0,1,n,N);

OpeL:=proc(a,b,c,n,N) local ope,i:
ope:=
-(b+n+1)*(a+n+1)*(a+b+2*n+4)+(2*n+3+a+b)*(a^2*c+a*b*c+2*a*c*n+2*b*c*n+2*c*n^2+a^2+2*a*b+3*a*c+4*a*n+b^2+3*b*c+4*b*n+6*c*n+4*n^2+6*a+6*b+4*c+12*n+8)*N-c^2*(n+2)*(2*n
+2+b+a)*(a+b+n+2)*N^2:

ope:=expand(ope/coeff(ope,N,2)):
add(factor(coeff(ope,N,i))*N^i,i=0..2):

end:





#AppxSeq1(ope,n,N,RE,K): It inputs a second order operator ope, in n where N is the shift operator N
#where it is known that ope(n,N) annihilates a sequence a(n) that goes to zero.
#It also inputs a pair RE=[[c0,c1],L] satisfying c0*a(0)+c1*a(1)=L where c0,c1 are integers and
#L is a rational number, It outputs the list of the first K terms in the rational approximations to a(0). Try
#AppxSeq1(N^2-N-1,n,N,[[11,12],3],20);
AppxSeq1:=proc(ope,n,N,RE,K) local c0,c1,L,gu,a0,lu,i,mu0,mu1:

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


c0:=RE[1][1]:
c1:=RE[1][2]:

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

L:=RE[2]:

gu:=SeqFromRec(ope,n,N,[a0,(L-c0*a0)/c1],K):

gu:=[op(3..nops(gu),gu)]:
lu:=[]:

for i from 2 to nops(gu) do
 mu0:=coeff(gu[i],a0,0):
 mu1:=coeff(gu[i],a0,1):
 if mu1=0 then
   RETURN(FAIL):
 fi:
 lu:=[op(lu),-mu0/mu1]:
od:
lu:
end:



#AppxSeq(a,b,c,K): The approximating sequence for  IntL(a,b,c,0); 
#Try:
#AppxSeq(0,0,1,40);
AppxSeq:=proc(a,b,c,K) local ope,n,N,RE:

RE:=FindRel(a,b,c):

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



ope:=OpeL(a,b,c,n,N):

AppxSeq1(ope,n,N,RE,K):

end:



#GLC1(a,b,c, K): A floating-point approximation to the constant
#IntGL(a,b,c,0) using the sequence to K terms 
#followed by the differene of two consecutive terms. Try:
#GLC1(0,0,1,100);
GLC1:=proc(a,b,c,K)  local RE, n,N,ope,lu:

RE:=FindRel(a,b,c):

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

ope:=OpeL(a,b,c,n,N):

lu:=AppxSeq1(ope,n,N,RE,K):

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


[evalf(lu[nops(lu)]),evalf(abs(lu[nops(lu)]-lu[nops(lu)-1]))]:

end:

#GLC(a,b,c): A floating-point approximation to the constant
#IntGL(a,b,c,0) supposed to be good for the number of digits
#Try:
#GLC(0,0,1);
GLC:=proc(a,b,c)  local lu,K,RE:


RE:=FindRel(a,b,c):

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

lu:=GLC1(a,b,c,100):

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

if lu[2]<10^(-3*Digits) then
 RETURN(lu[1]):
fi:

for K from 150 to 20000 by 50 do

lu:=GLC1(a,b,c,K):


if lu<>FAIL and lu[2]<10^(-3*Digits) then
 RETURN(lu[1]):
fi:

od:

FAIL:

end:


#IsGoodOpe(ope,n,N): Is ope(n,N) a good operator?
IsGoodOpe:=proc(ope,n,N) local gu:
gu:=denom(ope):
if subs(n=0,gu)<>0 then
 true:
else
 false:
fi:
end:

#AnBn(a,b,c,K): The first K terms in the sequences An and Bn in the expression
#IntGL(a,b,c,n)=B(n)*C-A(n)
#where C=IntGL(a,b,c,0); 
#It also returns the first K terms in the floating point of the sequence of errors
#Try:
#AnBn(0,0,1,50);
AnBn:=proc(a,b,c,K) local ope,n,N,RE,c0,c1,L,a0,i,gu,guF:

RE:=FindRel(a,b,c):
if RE=FAIL then
 RETURN(FAIL):
fi:


c0:=RE[1][1]:
c1:=RE[1][2]:

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

L:=RE[2]:

ope:=OpeL(a,b,c,n,N):

if not IsGoodOpe(ope,n,N)  then
 RETURN(FAIL):
fi:

gu:=SeqFromRec(ope,n,N,[a0,(L-c0*a0)/c1],K):

guF:=evalf(subs(a0=GLC(a,b,c),gu)):

[[seq(-coeff(gu[i],a0,0),i=1..nops(gu))],[seq(coeff(gu[i],a0,1),i=1..nops(gu))],guF]:

end:



#CnDn(a,b,c,K): The first K terms in the sequences gcd(numer(An[i]),number(Bn[i]))
#lcm(denom(An[i]),denom(Bn[i]))
#It also returns the last 20 terms of the logs of the normailizing sequence
#as well as the last 20 terms of the implied deltas
#Try:
#CnDn(0,0,1,100);
	CnDn:=proc(a,b,c,K) local gu,i,lu,vu,ku,beta1,beta2:

beta1:=Lcons(c)[1]:
beta2:=1/Lcons(c)[2]:

if K<20 then
 print(K, `shpild be at least 20 `):
 RETURN(FAIL):
fi:

gu:=AnBn(a,b,c,K):

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

lu:=
[
[seq(gcd(numer(gu[1][i]),numer(gu[2][i])),i=1..nops(gu[1]))],
[seq(lcm(denom(gu[1][i]),denom(gu[2][i])),i=1..nops(gu[1]))]
]:

vu:=evalf([seq((log(lu[2][i+1])-log(lu[1][i+1]))/i,i=K-20..K)]):

ku:=[ seq( (log(beta2)-vu[i])/(log(beta1)+vu[i]),i=1..nops(vu))]:

ku:=evalf(ku,10):
vu:=evalf(vu,10):

[lu,vu,ku]:

end:




#Lcons(c): The pair of numbers [beta1,beta2] such that An,Bn=Omega(beta1)^b, I(n)=A(n)-c*B(n)= Omega(beta2^n).
Lcons:=proc(c):
[(c+2+2*(c+1)^(1/2))/c^2, -(-c-2+2*(c+1)^(1/2))/c^2]
end:



#IR(L): IR(L): same as IntegerRelations[PSLQ](L)
IR:=proc(L): IntegerRelations[PSLQ](L):end:

rf1:=proc(a,n): (a+n-1)!/(a-1)!:end:

#IntGLsummand(a,b,c,n): The summand in the series representation for IntGL(a,b,c,0); Try:
#IntGLsummand(0,0,1,n);
IntGLsummand:=proc(a,b,c,n):

simplify((-1)^n*rf1(a+1,n)/rf1(a+b+2,n)*c^n):

end:


#IntGLsummandRF(a,b,c,n,RF): The summand in the series representation for IntGL(a,b,c) in terms of the raising-factorial
#denoted by RF. Try
#IntGLsummandRF(0,0,1,n,RF);
IntGLsummandRF:=proc(a,b,c,n,RF):

RF(a+1,n)/RF(a+b+2,n)*(-1)^n*c^n:

end:


#IntGL2F1(a,b,c): The summand in the series representation for IntGL(a,b,c,0) expressed as an
#2F1, a pair of lists [NumeratorParameters,DenominatorParameters], followed by the argument
#IntGL2F1(0,0,1);
IntGL2F1:=proc(a,b,c) local ku:
ku:=[[1,a+1],[a+b+2],-c]:
cat(`2F1(`,ku[1][1],`,`, ku[1][2], `;`, op(ku[2]), `;`, ku[3],`)`);
end:



###START PRIMES PROCEDURES

#Lc(n): the lcm of 1, ..., n. Try: Lc(1000);
Lc:=proc(n) local i:
lcm(seq(i,i=1..n)):
end:


#PpG(a,b,n,m1,m2,C,M): Given rational numbers a and b between 0 and 1 and positive integers m1,m2 and n
#and a subset C in {1,...M-1} and an integer M
#it is the product of p between m1 and m2  such that frac{n/p} is between a and b
#and p mod M is in C
PpG:=proc(a,b,n,m1,m2,C,M) local p,gu:
p:=2:

gu:=1:

p:=nextprime(m1):


while p<m2 do
 if (frac(n/p)>= a and frac(n/p)<b and member(p mod M, C)) then
 gu:=gu*p:
fi:
p:=nextprime(p):
od:

gu:

end:

#LcAB(n,A,B): the lcm of seq(A+B*i,i=1..n). Try: LcAB(1000,1,2);
LcAB:=proc(n,A,B) local i:
lcm(seq(A+B*i,i=1..n)):
end:


#PrimesF(N,n,k): inputs a HUGE integer N that is a product of small primes, and a much smaller integer n
#that generated N (via some process), outputs the list of length k, whose i-th entry is the list
#of primes>=sqrt(n) that show up with exponent i, followed by the list of primes that
#did not show up at all. Try:
#lu:=CnDn(1/3,1/3,1/15,1000)[1][2][-1]: PrimesF(lu,1000,6);

PrimesF:=proc(N,n,k) local gu,PR,KULAM,T,i,gadol,katan,j,p:

gu:=ifactors(N)[2]:

if nops(gu)=1 then
 RETURN(gu):
fi:

PR:={}:

for i from 1 to k do
T[i]:=[]:
od:

for j from 1 to nops(gu) do
 if gu[j][1]>=evalf(sqrt(n)) then
    T[gu[j][2]]:=[op(T[gu[j][2]]),gu[j][1]]:
  PR:=PR union {gu[j][1]}:
 fi:
od:

gadol:=max(op(PR)):
katan:=min(op(PR)):

KULAM:={}:

p:=nextprime(katan):
KULAM:={p}:

while p<=gadol do
 p:=nextprime(p):
 KULAM:=KULAM union {p}:
od:

PR:=KULAM minus PR:

[ [seq(T[i],i=1..k)],PR]:

end:



#Fs(S,n): Given a set S and a positive integer n finds sort([seq(evalf(frac(n/p),10),p in S)]);
Fs:=proc(S,n) local p:
sort([seq(evalf(frac(n/p),10),p in S)]);
end:


#Seg1(L,a,b): Given a  list L and integers a<b, extracts the segment between a and b. Try:
#Seg1([seq(i,i=1..100)],40,60);
Seg1:=proc(L,a,b) local L1,i,j:
L1:=sort(L):

if a>b then
 RETURN(FAIL):
fi:

for i from 1 while L1[i]<a do od:

for j from i while L[j]<b do od:

[op(i..j-1,L)]:
end:




#EvalPpG1(P,n1): evaluates a partial  PpG object at n=n1. 
#For GBC(0,0,0,0) (alias Zeta(2))
#EvalPpG1([[1, 1], [1, 2],[] ],1000  );

EvalPpG1:=proc(P,n1) local gu,lu,i,lu1:


if not (type(P,list) and nops(P)=3 and nops(P[1])=2 and nops(P[2])=2 and (P[3]=[] or {seq(nops(P[3][i]),i=1..nops(P[3]))}={6}) ) then
  print(`Bad input`):
  RETURN(FAIL):
fi:

gu:=P[1][1]^(P[1][2]*n1):

gu:=gu*Lc(P[2][1]*n1)^P[2][2]:
lu:=P[3]:

for i from 1 to nops(lu) do
lu1:=lu[i]:
gu:=gu*PpG(lu1[1],lu1[2],n1,trunc(lu1[3]*n1),trunc(lu1[4]*n1),lu1[5],lu1[6]    ):
od:

gu:

end:




#EvalPpG(PA,n1): evaluates a   PpG [Cn,Dn] where Cn and Dn are partial PpG objects
#For GBC(0,0,0,0) (alias Zeta(2))
#EvalPpG( [ [[1,1],[1,0],[]],    [[1, 1], [1, 2],[] ] ] ,1000  );
EvalPpG:=proc(PA,n1) :
if not (type(PA,list) and nops(PA)=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if EvalPpG1(PA[1],n1)=FAIL or EvalPpG1(PA[2],n1)=FAIL  then
 RETURN(FAIL):
fi:

 EvalPpG1(PA[2],n1)/EvalPpG1(PA[1],n1):
end:



#AsyPpG1(P): inputs a partial Pp expression  and outputs the
#EXACT value of limit of log( EvalPpG1(P,n))/n as n goes to infinity. Try
#AsyPpG1([[3, 3], [3, 1], [ [1/3,1,0,3/2,{0},1],  [1/3,1,3/2,3,{2},3],[2/3,1,0,3/2,{0},1]] ]  );
AsyPpG1:=proc(P) local gu,lu,i:


if not (type(P,list) and nops(P)=3 and nops(P[1])=2 and nops(P[2])=2 and 
(P[3]=[] or {seq(nops(P[3][i]),i=1..nops(P[3]))}={6} )) then
  print(`Bad input`):
  RETURN(FAIL):
fi:


gu:=P[1][2]*log(P[1][1]):


gu:=gu+P[2][1]*P[2][2]:


lu:=P[3]:

for i from 1 to nops(lu) do
gu:=gu+PpGlimit(op(lu[i])):
od:

gu:

end:

#AsyPpG(PA,n1): Inputs a GpP expression in the form [Cn,Dn] where Cn a, Dn are the partial PpG expressions for the sequences Cn and Dn (see paper)
#outputs nu sich that the INTEGERating sequence Dn/Cn is asympotic to exp(nu*n). 
#For GBC(0,0,0,0) (alias Zeta(2))
#AsyPpG(      [[1,1],[1,0],[]],    [[1, 1], [1, 2],[] ] ] ],1000  );
AsyPpG:=proc(PA) :
if not (type(PA,list) and nops(PA)=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if AsyPpG1(PA[1])=FAIL or AsyPpG1(PA[2])=FAIL  then
 RETURN(FAIL):
fi:

AsyPpG1(PA[2])-AsyPpG1(PA[1]):
end:


#The limit of log PpG(a,b,r1*n,r2*n,C,M)/n as n goes to infinity. Try
#PpGlimit(1/3,2/3,3/2,3,{2},3);
PpGlimit:=proc(a,b,r1,r2,C,M)
(PpLimit(a,b,r2)-PpLimit(a,b,r1))*nops(C)/numtheory[phi](M):
end:


#PrintPpG(P,n,L1,L2): P is the  same input as the P EvalPpG(P,n1) (q.v.) but n,L1 and L2 are symbolcs
#where n stands for n,  L1(n) stands for lcm(1...n) and L2(alpha,beta,a,b,C,M,n) stands for the
#product of all primes between a*n and b*n such that the fractional part of n/p is between alpha and beta
#and p mod M belongs to the set C where C is a subset of the set of nonneg. integers less than M relatively prime to it.
#If no modularity condition is present than M=1 and C={0}. Try:

#PrintPpG([[1, 1], [3, 1],[] ],n,LCM,P );
#PrintPpG([[3, 3], [3, 1], [ [1/3,1,0,3/2,{0},1],  [1/3,1,3/2,3,{2},3],[2/3,1,0,3/2,{0},1]] ],n,LCM,P );
#PrintPpG([[2, 6], [4, 1], [ [1/4,3/4,0,2,{0},1],  [3/4,1,0,4/3,{0},1], [3/4,1,0,4/3,{0},1] ] ],n, LCM,P );
PrintPpG:=proc(P,n,L1,L2) local gu,lu,i,lu1:


if not (type(P,list) and nops(P)=3 and nops(P[1])=2 and nops(P[2])=2 and (P[3]=[] or {seq(nops(P[3][i]),i=1..nops(P[3]))}={6}) ) then
  print(`Bad input`):
  RETURN(FAIL):
fi:

gu:=1/P[1][1]^(P[1][2]*n):

gu:=gu*L1(P[2][1]*n)^P[2][2]:
lu:=P[3]:

for i from 1 to nops(lu) do
lu1:=lu[i]:
gu:=gu*L2(op(lu1),n):
od:

gu:

end:


####################
# PpLimit(a, b, r): Return the limit of log(Pp(a, b, n, rn)) / n as n tends to
# infinit.
# Try:
#       evalf([log(Pp(1/4, 3/4, 1000, 2 * 1000)) / 1000, PpLimit(1/4, 3/4, 2)]);
#       evalf([log(Pp(3/4, 1, 1000, 4 * 1000 / 3)) / 1000, PpLimit(3/4, 1, 4/3)]);
#       evalf([log(Pp(1/10, 1/2, 1000, 5 * 1000)) / 1000, PpLimit(1/10, 1/2, 5)]);
#       evalf([log(Pp(1/10, 2/10, 2000, 3 * 2000)) / 2000, PpLimit(1/10, 2/10, 3)]);
PpLimit:=proc(a, b, r)
if r=0 then
 RETURN(0):

elif a=0 and b=1 then
r:
elif a=0 then
r-PpLimit(b, 1, r):

else
    Psi(b) - Psi(a) - ifelse(a < 1 / r, 1 / a - max(r, 1 / b), 0);
fi:
end:



#PrintPpG1(P,n,L1,L2): P is the  same input as the P EvalPpG1(P,n1) (q.v.) but n,L1 and L2 are symbolcs
#where n stands for n,  L1(n) stands for lcm(1...n) and L2(alpha,beta,a,b,C,M,n) stands for the
#product of all primes between a*n and b*n such that the fractional part of n/p is between alpha and beta
#and p mod M belongs to the set C where C is a subset of the set of nonneg. integers less than M relatively prime to it.
#If no modularity condition is present than M=1 and C={0}. Try:
PrintPpG1:=proc(P,n,L1,L2) local gu,lu,i,lu1:


if not (type(P,list) and nops(P)=3 and nops(P[1])=2 and nops(P[2])=2 and (P[3]=[] or {seq(nops(P[3][i]),i=1..nops(P[3]))}={6}) ) then
  print(`Bad input`):
  RETURN(FAIL):
fi:

gu:=P[1][1]^(P[1][2]*n):

gu:=gu*L1(P[2][1]*n)^P[2][2]:
lu:=P[3]:

for i from 1 to nops(lu) do
lu1:=lu[i]:
gu:=gu*L2(op(lu1),n):
od:

gu:

end:


#PrintPpG(PA,n,L1,L2): P is the  same input as the P EvalPpG(P,n1) (q.v.) but n,L1 and L2 are symbolcs
#where n stands for n,  L1(n) stands for lcm(1...n) and L2(alpha,beta,a,b,C,M,n) stands for the
#product of all primes between a*n and b*n such that the fractional part of n/p is between alpha and beta
#and p mod M belongs to the set C where C is a subset of the set of nonneg. integers less than M relatively prime to it.
#If no modularity condition is present than M=1 and C={0}. Try:
PrintPpG:=proc(PA,n,L1,L2) :

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

if PrintPpG1(PA[1],n,L1,L2)=FAIL or PrintPpG1(PA[2],n,L1,L2)=FAIL then
 RETURN(FAIL):
fi:

PrintPpG1(PA[2],n,L1,L2)/PrintPpG1(PA[1],n,L1,L2):
end:


###END PRIMES PROCEDURES

#MyIDs(C,F,F0,N): Given a constant C in decimals and another constant
#let's call it F, whose floating-point if F0
# and a positive integer N
#tries to express C as (a*F+b)/(c*F+d) for a,b,c,d from -N to N using PSLQ
#MyIDs(evalf((log(2)-2)/(2*log(2)+3)),log(2),evalf(log(2)),100);
MyIDs:=proc(C,F,F0,N) local gu,mu,i:

gu:=IntegerRelations[PSLQ]([1,evalf(C),evalf(F0),evalf(C*F0)]):


if max(seq(abs(gu[i]),i=1..4)) >N then
 RETURN(FAIL):
fi:

mu:=-(gu[1]+gu[3]*F)/(gu[2]+gu[4]*F):

if abs(evalf(subs(F=F0,mu)-C))>1/10^(Digits-7) then
  RETURN(FAIL):
fi:

mu:

end:


#TheoremL(a,b,c,K,P): Inputs rational numbers a,b,c, and a large positive integer K and outputs
#a theorem regarding the constant
#print(Int(x^(a)*(1-x)^(b)/(1+c*x),x=0..1)/Int(x^(a)*(1-x)^(b),x=0..1).
#Either a suggested proof of irrationality or a way to compute it exponentially  fast.
#P is 0 if no INTEGERating factor is known, otherwise it is a proposed one.
#Try:
#TheoremL(0,0,1,2000,[ [1,0],[1,1],[]] ):

TheoremL:=proc(a,b,c,K,P) local n,N,gu,x,y,ope,lu,X,A,B,beta1,beta2,F,C,E,A1,B1,delta,d1,ka,nuE,deltE,WADIM,su2F1,deltE1:
beta1:=Lcons(c)[1]:
beta2:=1/Lcons(c)[2]:

gu:=CnDn(a,b,c,K):
lu:=AnBn(a,b,c,K):
ka:=GCL(a,b,c):
su2F1:=IntGL2F1(a,b,c,n):

 ope:=OpeL(a,b,c,n,N):

print(``):

if min(op(gu[3]))<0 then
 print(`Very fast Computation of rational approximations to the constant `):
 print( su2F1 ):
 print(``):
 print(`By Shalosh B. Ekhad `):
 print(``):


 print(`Theorem: , let A(n), B(n), be two sequences of rational numbers that satisfy the second-order recurrence`):
 print(``):
 print( add(coeff(ope,N,i)*X(n+i),i=0..degree(ope,N)) =0):
 print(``):
 print(`Subject to the initial conditions`):
  print(``):
 print(A(0)=lu[1][1],A(1)=lu[1][2]):
  print(``):
  print(B(0)=lu[2][1],B(1)=lu[2][2]):
 print(``):
 print(`Then`, A(n)/B(n), `approximates the constant of the title, let's call it C`, cat(`2F1`,`(`, su2F1,`)` ) ):
print(`That is easily seen to be given by the following Alladi-Robinson-type  integral`):
 print(``):
  print( C=Int(x^(a)*(1-x)^(b)/(1+c*x),x=0..1)/Beta(1+a,1+b)):


 print(``):
  print(`with an error that is OMEGA of`, 1/(beta1*beta2)^n, ` that in floating point is`,  1/(beta1*beta2)^n   ):
 print(``):


if identify(ka)<>ka then
  print(`Comment: Note that this constant appears to be `, identify(ka) ):
  print(`Prove it!`):
fi:


if P=0 then

print(`Comment: while this sequence does not lead to an irrationality proof, for the record, the delta, happens to be roughly`):
print(``):
print(min(op(gu[3]))):

print(``):

else
nuE:=AsyPpG(P):


deltE1:=  (log(beta2)-nuE)/(log(beta1)+nuE):
print(`Comment: while this sequence does not lead to an irrationality proof, for the record, the delta, happens to be exactly`):
print(``):
print(deltE1):
print(``):
print(`that in decimals is`):

print(evalf(deltE1,10)):


print(`In order to prove the exact value of the delta such that there exist integer sequences An, Bn such that abs(c-An/Bn)<CONST/Bn^(1+delta)`):
print(`we need the following lemma`):


print(`Lemma : Let E(n) be`):
print(``):
print(PrintPpG(P,n,LCM,PP)):
print(``):
if P[1][2][2]<>0 or P[2][2][2]<>0 then
print(`where LCM(1..m) is the least-common-multiple of the first m integers, and`):
fi:
if P[1][3]<>[] or P[2][3]<>[] then
print(`where for 0<alpha<beta<1, a,b, an integer M and a subset C of {0,...M-1}`):
print(` PP(alpha,beta,a,b,C,M,n) is the product of all r primes between a*n and b*n, such that of you take them mod M it belongs to C`):
print(`and the fractional part of n/p is between alpha and beta`):
fi:
print(``):
print(`Then there exists a sequence G(n) of integers that is asymptotically negligible, i.e. the limit of log( G(n) )/n tends to 0 as n goes to infinity`):
print(`such that `):
 print(` A1(n):=E(n)G(n)A(n), B1(n):=E(n)G(n)B(n) are BOTH integers`):
 print(``):
print(`For example, G(n) for n=`, K, `equals ` ):
WADIM:=lcm(denom(EvalPpG(P,K)*lu[1][-1]),denom(EvalPpG(P,K)*lu[2][-1])):
print(ifactor(WADIM)):
print(`whose log divided by n is`):
print(``):
print(evalf(log(WADIM)/K,10)):
print(``):
print(`Indeed peanuts. `):
print(``):
 print(`Using standard asymptotics derived from the Prime Number Theorem it is readily seen that the limit of`):
print(` log(E(n))/n , and hence log(E(n)*G(n))/n, equals exactly`, nuE ):
 print(``):
print(`that equals in decimals`, evalf(nuE,10)):
print(``):  
 print(`Just as a sanity check, the empircal numerical values of nu for E(n) from`, K+1-nops(gu[2]), `to `, K+1, `are `):
 print(``):
 print(gu[2]):
print(``):
 print(`Multiplying F(n) by E(n)*G(n) we get `):
print(``):
print( E(n)*F(n)*G(n)=B1(n)*c-A1(n)):
print(``):
print(`and this implies that `):
print(``):
print(abs(c-A1(n)/B1(n)) <= CONSTANT/B1(n)^(1+delta)):
print(``):

print(` where `, delta= (log(beta2)-nuE)/(log(beta1) + nuE) ):
print(``):
print(`that in decimals is, `,  evalf( (log(beta)-nuE)/(log(beta) + nuE) ,10)   )  :
print(``):
print(`Unfortunaley, delta is  negative `):
print(``):
print(`We leave it to the reader to fill-in the details.`):
print(``):

fi:
	

 print(`Proof, consider the Alladi-Robinson type-integral `):
 print(F(n)=Int(x^a*(1-x)^b/(1+c*x)*(x*(1-x)/(1+c*x))^n,x=0..1)/Beta(a+1,b+1) ):
 print(``):
 print(`Then `, F(0)=B(0)*c-A(0), F(1)=c*B(1)-A(1) ):
 print(``):
 print(`and F(n) also satisfies the above recurrence, thanks to the amazing classical Almkvist-Zeilberger algorithm `):
 print(``):
  print(`Hence`, F(n)=B(n)*c-A(n)):
 print(``):
  print(`By a simple bound of the integrand,  F(n) is OMEGA of`, 1/beta2^n, ` and  by the Poincare lemma, B(n) (and for that matter, A(n)) are OMEGA of`, beta1^n):
 print(``):
 print(`Dividing by B(n) gives that A(n)/B(n)-c is OMEGA of `, 1/(beta1*beta2)^n , `QED. `):



else


 print(`Sketch of an Irrationality Proof of the constant `):
print(su2F1 ):

  print(``):
 print(`By Shalosh B. Ekhad `):
   print(``):
 print(`Theorem: The constant of the title `):
   print(``):
print(su2F1 ):
print(``):
print(`let's call it C, that is easily seen to be given by the following Alladi-Robinson-type integral`):
print(``):
 print( C=Int(x^(a)*(1-x)^(b)/(1+c*x),x=0..1)/Beta(1+a,1+b) ):
    print(``):

if P=0 then
 print(`is irrational, with an irrationality measure`, 1+ (log(beta1)+nu)/(log(beta2)-nu), `for a certain number nu `):
print(`that is approximately `, max(op(gu[2])) , ` yielding an irrationality measure that is approximately `, evalf(1+1/min(op(gu[3])),10)):
print(``):

 print(`We hope that the reader can find  nu exactly. `):

else

nuE:=AsyPpG(P):

deltE:= 1+ (log(beta1)+nuE)/(log(beta2)-nuE):

 print(`is irrational, with an irrationality measure`, deltE):
print(`that in decimals equals`, evalf( deltE,10) ):
print(``):
fi:

if identify(ka)<>ka then
  print(`Comment: Note that this constant appears to be `, identify(ka) ):
  print(`Prove it!`):
fi:

 print(`We need two  lemmas `):
    print(``):

 print(`Lemma ONE: , let A(n), B(n), be two sequences of rational numbers that satisfy the second-order recurrence`):
 print(``):
 print( add(coeff(ope,N,i)*X(n+i),i=0..degree(ope,N)) =0):
 print(``):
 print(`Subject to the initial conditions`):
  print(``):
 print(A(0)=lu[1][1],A(1)=lu[1][2]):
  print(``):
  print(B(0)=lu[2][1],B(1)=lu[2][2]):
 print(``):
 print(`Then`, A(n)/B(n), `approximates the constant of the title`, cat(`2F1`,`(`, su2F1, `)` )):
print(` let's call it C, that can be easily seen to be given by the Alladi-Robinson-type integral`):
 print(``):
  print( C=Int(x^(a)*(1-x)^(b)/(1+c*x),x=0..1)/Beta(1+a,1+b) ):
 print(``):
  print(`with an error that is OMEGA of`, (1/beta1*beta2)^n, ` that in floating point is`,  evalf((1/(beta1*beta2))^n) ):
 print(``):
 print(`Proof, consider the Alladi-Robinson type-integral `):
 print(F(n)=Int(x^a*(1-x)^b/(1+c*x)*(x*(1-x)/(1+c*x))^n,x=0..1)/Beta(1+a,1+b)):
 print(``):
 print(`Then `, F(0)=B(0)*c-A(0), F(1)=c*B(1)-A(1) ):
 print(``):
 print(`and F(n) also satisfies the above recurrence, thanks to the amazing multivariable Almkvist-Zeilberger algorithm `):
 print(``):
  print(`Hence`, F(n)=B(n)*c-A(n)):
 print(``):
  print(`By a simple bound of the integrand,  F(n) is OMEGA of`, 1/beta2^n, ` and  by the Poincare lemma, B(n) (and for that matter, A(n)) are OMEGA of`, beta1^n):
 print(``):
 print(`Dividing by B(n) gives that A(n)/B(n)-c is OMEGA of `, 1/(beta1*beta2)^n , `QED. `):
print(``):

 print(`we now claim that the sequence of RATIONAL numbers A(n),B(n), can be multiplied by another sequence of rational numbers `):
 print(`E(n) such that both A(n)E(n) and B(n)E(n) are integers `):
 print(``):


if P=0 then
 print(`Lemma TWO: There exists a sequence of rational numbers, whose prime factorizations consists of small primes, that hopefully`):
 print(`can be described (and proved) explicity, that we leave to the expert reader such that `):
 print(` A1(n):=E(n)A(n), B1(n):=E(n)B(n) are BOTH integers`):
 print(``):
 print(`Furthermore there exists a contant, nu, that hopefully the learned reader can determine such that E(n) is OMEGA of `, exp(nu*n) ):
 print(``):
 print(`The empircal values of nu for E(n) from`, K+1-nops(gu[2]), `to `, K+1, `are `):
 print(``):
 print(gu[2]):
print(``):
 print(`Multiplying F(n) by E(n) we get `):
print(``):
print( E(n)*F(n)=B1(n)*c-A1(n)):
print(``):
print(`and this implies that `):
print(``):
print(abs(c-A1(n)/B1(n)) <= CONSTANT/B1(n)^(1+delta)):
print(``):
print(` where `, delta= (log(beta2)-nu)/(log(beta1) + nu) ):
print(``):
print(`Using the above values of nu for E(n) from`, K+1-nops(gu[3]), `to `, K+1, `the estimated deltas are `):
print(``):
print(gu[3]):
print(``):
print(`As you can see, they are all positive `):
print(``):
print(`We leave it to the reader to fill-in the details.`):
print(``):


else

print(`Lemma TWO: Let E(n) be`):
print(``):
print(PrintPpG(P,n,LCM,PP)):
print(``):
if P[1][2][2]<>0 or P[2][2][2]<>0 then
print(`where LCM(1..m) is the least-common-multiple of the first m integers, and`):
fi:
if P[1][3]<>[] or P[2][3]<>[] then
print(`where for 0<alpha<beta<1, a,b, an integer M and a subset C of {0,...M-1}`):
print(` PP(alpha,beta,a,b,C,M,n) is the product of all r primes between a*n and b*n, such that of you take them mod M it belongs to C`):
print(`and the fractional part of n/p is between alpha and beta`):
fi:
print(``):
print(`Then there exists a sequence G(n) of integers that is asymptotically negligible, i.e. the limit of log( G(n) )/n tends to 0 as n goes to infinity`):
print(`such that `):
 print(` A1(n):=E(n)G(n)A(n), B1(n):=E(n)G(n)B(n) are BOTH integers`):
 print(``):
print(`For example, G(n) for n=`, K, `equals ` ):
WADIM:=lcm(denom(EvalPpG(P,K)*lu[1][-1]),denom(EvalPpG(P,K)*lu[2][-1])):
print(ifactor(WADIM)):
print(`whose log divided by n is`):
print(``):
print(evalf(log(WADIM)/K,10)):
print(``):
print(`Indeed peanuts. `):
print(``):
 print(`Using standard asymptotics derived from the Prime Number Theorem it is readily seen that the limit of`):
print(` log(E(n))/n , and hence log(E(n)*G(n))/n, equals exactly`, nuE ):
 print(``):
print(`that equals in decimals`, evalf(nuE,10)):
print(``):  
 print(`Just as a sanity check, the empircal numerical values of nu for E(n) from`, K+1-nops(gu[2]), `to `, K+1, `are `):
 print(``):
 print(gu[2]):
print(``):
 print(`Multiplying F(n) by E(n)*G(n) we get `):
print(``):
print( E(n)*F(n)*G(n)=B1(n)*c-A1(n)):
print(``):
print(`and this implies that `):
print(``):
print(abs(c-A1(n)/B1(n)) <= CONSTANT/B1(n)^(1+delta)):
print(``):


print(` where `, delta= (log(beta2)-nuE)/(log(beta1) + nuE) ):
print(``):
print(`that in decimals is, `,  evalf( (log(beta2)-nuE)/(log(beta1) + nuE) ,10)   )  :
print(``):
print(`As you can see, delta is  positive `):
print(``):
print(`We leave it to the reader to fill-in the details.`):
print(``):



fi:
fi:

end:


# Search(N, K): All the triples [a,b,c] = [a'/N, b'/N, c'/N] with a' and b'
# between 0 and N-1, and c' between 1 and N such that the deltas for
# CnDn(a,b,c) seems to be positive.
# Try Search(2,30);
Search:=proc(N,K)
    option remember:
    local a, b, c, params, seen, good, cndn, deltas,mu:
    seen:={}:
    good := []:
    for a from 0 to N-1 do
    for b from 0 to N-1 do
    for c from 1 to N do
        params := [a / N, b / N, c / N]:

        if member(params, seen) then
            next:
        fi:

        seen := seen union {params}:

        if lcm(op(denom(params))) <> N then
            next:
        fi:

        cndn := CnDn(op(params), K):

        if mu = FAIL then
            next:
        fi:

        deltas := cndn[3]:
        if min(deltas) > 0 then
            good := [op(good), params]:
        fi:
    od:
    od:
    od:

    good:
end:

Hopefuls1:=[[[0, 0, 1]], [[0, 0, 1/2], [1/2, 1/2, 1/2], [1/2, 1/2, 1]], [[0, 0, 1/3], [0, 0, 2/3], [1/3, 1/3, 1/3
], [1/3, 1/3, 1], [1/3, 2/3, 1/3], [1/3, 2/3, 2/3], [1/3, 2/3, 1], [2/3, 1/3, 1/3], [2/3, 1/3, 2/3], [
2/3, 1/3, 1], [2/3, 2/3, 1/3], [2/3, 2/3, 1]], [[0, 0, 1/4], [0, 1/2, 1/4], [1/4, 1/4, 1/4], [1/4, 1/4
, 1/2], [1/4, 1/4, 1], [1/4, 3/4, 1/4], [1/4, 3/4, 1/2], [1/4, 3/4, 1], [1/2, 1/4, 1/4], [1/2, 1/2, 1/
4], [1/2, 1/2, 3/4], [1/2, 3/4, 1/4], [3/4, 1/4, 1/4], [3/4, 1/4, 1/2], [3/4, 1/4, 1], [3/4, 3/4, 1/4]
, [3/4, 3/4, 1/2], [3/4, 3/4, 1]], [[0, 0, 1/5], [0, 0, 2/5], [0, 0, 3/5], [0, 1/5, 1/5], [0, 2/5, 1/5
], [0, 3/5, 1/5], [0, 4/5, 1/5], [1/5, 1/5, 1/5], [1/5, 4/5, 1/5], [1/5, 4/5, 2/5], [1/5, 4/5, 1], [2/
5, 2/5, 1/5], [2/5, 3/5, 1/5], [2/5, 3/5, 2/5], [2/5, 3/5, 1], [3/5, 2/5, 1/5], [3/5, 2/5, 2/5], [3/5,
2/5, 1], [3/5, 3/5, 1/5], [4/5, 1/5, 1/5], [4/5, 1/5, 2/5], [4/5, 1/5, 1], [4/5, 4/5, 1/5]], [[0, 0, 1
/6], [0, 1/3, 1/6], [0, 1/2, 1/6], [0, 2/3, 1/6], [1/6, 1/6, 1/6], [1/6, 1/6, 1/3], [1/6, 1/6, 1/2], [
1/6, 1/6, 2/3], [1/6, 1/6, 1], [1/6, 1/3, 1/6], [1/6, 1/3, 1/3], [1/6, 5/6, 1/6], [1/3, 1/6, 1/6], [1/
3, 1/3, 1/6], [1/3, 1/3, 1/2], [1/3, 2/3, 1/6], [1/3, 2/3, 1/2], [1/2, 0, 1/3], [1/2, 1/6, 1/6], [1/2,
1/2, 1/6], [1/2, 1/2, 1/3], [1/2, 1/2, 2/3], [1/2, 1/2, 5/6], [1/2, 5/6, 1/6], [2/3, 1/3, 1/6], [2/3,
1/3, 1/2], [2/3, 2/3, 1/6], [2/3, 2/3, 1/2], [2/3, 5/6, 1/6], [5/6, 1/6, 1/6], [5/6, 2/3, 1/3], [5/6,
5/6, 1/6], [5/6, 5/6, 1/3], [5/6, 5/6, 1/2], [5/6, 5/6, 2/3], [5/6, 5/6, 1]], [[0, 0, 1/7], [0, 0, 2/7
], [0, 0, 3/7], [0, 1/7, 1/7], [0, 2/7, 1/7], [0, 3/7, 1/7], [0, 4/7, 1/7], [0, 5/7, 1/7], [0, 6/7, 1/
7], [1/7, 1/7, 1/7], [1/7, 6/7, 1/7], [1/7, 6/7, 2/7], [1/7, 6/7, 1], [2/7, 2/7, 1/7], [2/7, 5/7, 1/7]
, [2/7, 5/7, 2/7], [2/7, 5/7, 1], [3/7, 3/7, 1/7], [3/7, 4/7, 1/7], [3/7, 4/7, 2/7], [3/7, 4/7, 1], [4
/7, 3/7, 1/7], [4/7, 3/7, 2/7], [4/7, 3/7, 1], [4/7, 4/7, 1/7], [5/7, 2/7, 1/7], [5/7, 2/7, 2/7], [5/7
, 2/7, 1], [5/7, 5/7, 1/7], [6/7, 1/7, 1/7], [6/7, 1/7, 2/7], [6/7, 1/7, 1], [6/7, 6/7, 1/7]], [[0, 0,
1/8], [0, 0, 3/8], [0, 1/4, 1/8], [0, 1/2, 1/8], [0, 3/4, 1/8], [1/8, 1/8, 1/8], [1/8, 1/8, 1/4], [1/8
, 3/8, 1/8], [1/8, 7/8, 1/8], [1/8, 7/8, 1/4], [1/4, 1/4, 1/8], [1/4, 1/2, 1/8], [1/4, 3/4, 1/8], [3/8
, 1/8, 1/8], [3/8, 3/8, 1/8], [3/8, 3/8, 1/4], [3/8, 5/8, 1/8], [3/8, 5/8, 1/4], [1/2, 0, 1/8], [1/2,
1/8, 1/8], [1/2, 1/4, 1/8], [1/2, 3/8, 1/8], [1/2, 1/2, 1/8], [1/2, 1/2, 3/8], [1/2, 1/2, 5/8], [1/2,
5/8, 1/8], [1/2, 3/4, 1/8], [1/2, 7/8, 1/8], [5/8, 3/8, 1/8], [5/8, 3/8, 1/4], [5/8, 5/8, 1/8], [5/8,
5/8, 1/4], [5/8, 7/8, 1/8], [3/4, 1/4, 1/8], [3/4, 1/2, 1/8], [3/4, 3/4, 1/8], [7/8, 1/8, 1/8], [7/8,
1/8, 1/4], [7/8, 5/8, 1/8], [7/8, 7/8, 1/8], [7/8, 7/8, 1/4]], [[0, 0, 1/9], [0, 0, 2/9], [0, 0, 4/9],
[0, 1/3, 1/9], [0, 1/3, 2/9], [0, 2/3, 1/9], [0, 2/3, 2/9], [1/9, 1/9, 1/9], [1/9, 2/9, 1/9], [1/9, 5/
9, 1/9], [1/9, 8/9, 1/9], [1/9, 8/9, 2/9], [1/9, 8/9, 1/3], [2/9, 1/9, 1/9], [2/9, 2/9, 1/9], [2/9, 4/
9, 1/9], [2/9, 7/9, 1/9], [2/9, 7/9, 2/9], [2/9, 7/9, 1/3], [1/3, 1/9, 1/9], [1/3, 2/9, 1/9], [1/3, 1/
3, 1/9], [1/3, 1/3, 2/9], [1/3, 4/9, 1/9], [1/3, 5/9, 1/9], [1/3, 2/3, 1/9], [1/3, 2/3, 2/9], [1/3, 7/
9, 1/9], [1/3, 8/9, 1/9], [4/9, 2/9, 1/9], [4/9, 4/9, 1/9], [4/9, 5/9, 1/9], [4/9, 5/9, 2/9], [4/9, 5/
9, 1/3], [4/9, 8/9, 1/9], [5/9, 1/9, 1/9], [5/9, 4/9, 1/9], [5/9, 4/9, 2/9], [5/9, 4/9, 1/3], [5/9, 5/
9, 1/9], [5/9, 7/9, 1/9], [2/3, 1/9, 1/9], [2/3, 2/9, 1/9], [2/3, 1/3, 1/9], [2/3, 1/3, 2/9], [2/3, 1/
3, 4/9], [2/3, 4/9, 1/9], [2/3, 5/9, 1/9], [2/3, 2/3, 1/9], [2/3, 2/3, 2/9], [2/3, 8/9, 1/9], [7/9, 2/
9, 1/9], [7/9, 2/9, 2/9], [7/9, 2/9, 1/3], [7/9, 5/9, 1/9], [7/9, 7/9, 1/9], [7/9, 8/9, 1/9], [8/9, 1/
9, 1/9], [8/9, 1/9, 2/9], [8/9, 1/9, 1/3], [8/9, 4/9, 1/9], [8/9, 7/9, 1/9], [8/9, 8/9, 1/9]]]:



#Info1(a,b,c,K): inputs a triple (a,b,c) and outputs a list consisting of
#(i) the constant as a 2F1, (ii) the floating point to 100 digits,
#(iii) the recurrence operator annihilating An and Bn (iv) the pair of initial conditions for An and Bn
#i.e. [[A(0),A(1)],[B(0),B(1)]] (v) the estimated value of the limit of the log of the INTEGERATING functor divided by n (vi) the
#estimated implied irrationality measure using K values (vii) an identification if successful. Try:
#Info1(0,0,1,2000);
Info1:=proc(a,b,c,K) local mu,lu,n,N:

mu:=AnBn(a,b,c,2):
lu:=CnDn(a,b,c,K):
[IntGL2F1(a,b,c),evalf(GLC(a,b,c),100),OpeL(a,b,c,n,N),[[op(1..2,mu[1])],[op(1..2,mu[2])]],lu[2][-1],evalf(1+1/lu[3][-1],10),[a,b,c]]:
end:


#Hopefuls1v(K): Lots of constants that probably have an irrationality proof, together for each the relevant information. Try:
#Hopefuls1v(2000):
Hopefuls1v:=proc(K) local gu,i,j,L,ka,ka1,lu:
gu:=Hopefuls1():

print(`This is a collection of candidates for Apery-style proofs using generalized Alladi-Robinson integrals`):
print(``):

print(`Each entry is a list of length 7 where `):
print(``):
print(`The first entry is the definition of the constant as a 2F1`):
print(``):
print(`The second entry is its value to 100 decimals`):
print(``):
print(`The third entry is the second-order linear recurrence operator in n and N  annihilating Both sequence of rational numbers A(n) and B(n)`):
print(`such that that A(n)/B(n) converges to the constant`):
print(``):
print(`The fourth entry is the pair [[A(0),A[1]],[B(0),B(1)]], of initial conditions enabling the fast computations of A(n) and B(n) using the recurrence`):
print(``):
print(`The fifth entry is the estimate for the limit of log(E(n))/n as n goes to infinity of the INTEGERating factor, i.e. the sequence`):
print(`of rational numbers such that A(n)*E(n) and B(n)*E(n) are BOTH integers`):
print(``):
print(`The sixth entry is the estimate for the implied irrationality measure.`):
print(``):
print(``):
print(`The seventh entry is the triple [a,b,c] such that the constant is Int(x^a*(1-x)^b/(1+c*x),x=0..1)/Beta(a+1,b+1).`):
print(``):

print(`If Maple can identify the constant, then there is a comment to that effect.`):
print(``):



for i from 1 to nops(gu) do
print(``):
print(`denominator `, i):
L:=gu[i]:

for j from 1 to nops(L) do
 lu:=L[j]:

if min(op(CnDn(op(lu),K)[3]))>0 then
 lprint(Info1(op(lu),K)):

 ka:=GLC(op(lu)):
 ka1:=identify(ka):
  if ka<>ka1 then  
  print(`Comment: This probably equals`, ka1):
 fi:
fi:

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

end:


####start Generalized procedures


#GLC1g(a,b,c, d,K): A floating-point approximation to the constant
#IntGLg(a,b,c,d,0) using the sequence to K terms 
#followed by the differene of two consecutive terms. Try:
#GLC1g(0,0,1,0,100);
GLC1g:=proc(a,b,c,d,K)  local RE, n,N,ope,lu:

RE:=FindRelG(a,b,c,d):

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

ope:=OpeLg(a,b,c,d,n,N):

lu:=AppxSeq1(ope,n,N,RE,K):

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


[evalf(lu[nops(lu)]),evalf(abs(lu[nops(lu)]-lu[nops(lu)-1]))]:

end:

#GLCg(a,b,c,d): A floating-point approximation to the constant
#IntGLg(a,b,c,d,0) supposed to be good for the number of digits
#Try:
#GLCg(0,0,1,0);
GLCg:=proc(a,b,c,d)  local lu,K,RE:


RE:=FindRelG(a,b,c,d):

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

lu:=GLC1g(a,b,c,d,100):

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

if lu[2]<10^(-3*Digits) then
 RETURN(lu[1]):
fi:

for K from 150 to 20000 by 50 do

lu:=GLC1g(a,b,c,d,K):


if lu<>FAIL and lu[2]<10^(-3*Digits) then
 RETURN(lu[1]):
fi:

od:

FAIL:

end:



#OpeLg:=proc(a,b,c,d,n,N): The operator in n and the shift operator N, annihaling int(x^a*(1-x)^b/(1+c*x)^d (x*(1-x)/(1+c*d))^n,x=0..1):
#Try:
#OpeLg(0,0,1,0,n,N);
OpeLg:=proc(a,b,c,d,n,N) local ope:
ope:=
-(b+n+1)*(a+n+1)*(a+b+2*n+4)+(2*n+3+a+b)*(a^2*c+a*b*c-a*c*d+2*a*c*n+b*c*d+2*b*c*n+2*c*n^2+a^2+2*a*b+3*a*c+4*a*n+b^2+3*b*c+4*b*n+6*c*n+4*n^2+6*a+6*b+4*c+12*n+8)*N-c^
2*(d+n+2)*(2*n+2+b+a)*(a+b-d+n+2)*N^2:

ope:=ope/coeff(ope,N,2):
add(factor(coeff(ope,N,i))*N^i,i=0..2):
end:


#IntGLg(a,b,c,d,n1);
#IntGLg(0,0,1,0,0);
IntGLg:=proc(a,b,c,d,n) local x,y,lu,mu:
lu:=int(x^a*(1-x)^b*x^n*(1-x)^n/(1+c*x)^(n+d+1),x=0..1):
mu:=int(x^a*(1-x)^b/(1+c*x)^d,x=0..1):
lu/mu:
end:


#FindRelG(a,b,c,d): A  relationship
#between IntGLg(a,b,c,d,0) and  IntGLg(a,b,c,d,1) .  The output is of the form [[a0,a1],c] , that means that
#a0*IntGLg(a,b,c,d,0)+a1*IntGLg(a,b,c,d,1)=c. Try:
#FindRelG(0,0,1,0);
FindRelG:=proc(a,b,c,d) local x,c0,c1,gu,eq,var,lu,A,B,C,gu1,gu2,hal,ka,i:

option remember:

lu:=x^(a)*(1-x)^(b)/(1+c*x)^d:

gu1:=normal(diff((c0+c1*x)*x*(1-x)*lu/(1+c*x),x)):

gu1:=normal(gu1/lu):

gu2:=normal(A/(1+c*x) +B*x*(1-x)/(1+c*x)^2-C):

gu:=normal(gu1-gu2):

gu:=numer(gu):

var:={c0,c1,A,B,C}:

eq:={seq(coeff(gu,x,i),i=0..degree(gu,x))}:

var:=solve(eq,var):

hal:=subs(var,[A,B,C]):

if hal[1]=0 then
 RETURN(FAIL):
fi:

hal:=[1,normal(hal[2]/hal[1]),normal(hal[3]/hal[1])]:

hal:

ka:=lcm(seq(denom(hal[i]),i=1..nops(hal))):

hal:=ka*hal:

[[hal[1],hal[2]],hal[3]]:

end:


#AnBnG(a,b,c,d,K): The first K terms in the sequences An and Bn in the expression
#IntGLg(a,b,c,d,n)=B(n)*C-A(n)
#where C=IntGLg(a,b,c,d,0); 
#It also returns the first K terms in the floating point of the sequence of errors
#Try:
#AnBnG(0,0,1,0,50);
AnBnG:=proc(a,b,c,d,K) local ope,n,N,RE,c0,c1,L,a0,i,gu,guF:

RE:=FindRelG(a,b,c,d):
if RE=FAIL then
 RETURN(FAIL):
fi:


c0:=RE[1][1]:
c1:=RE[1][2]:

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

L:=RE[2]:

ope:=OpeLg(a,b,c,d,n,N):

if not IsGoodOpe(ope,n,N)  then
 RETURN(FAIL):
fi:

gu:=SeqFromRec(ope,n,N,[a0,(L-c0*a0)/c1],K):

guF:=evalf(subs(a0=GLCg(a,b,c,d),gu)):

[[seq(-coeff(gu[i],a0,0),i=1..nops(gu))],[seq(coeff(gu[i],a0,1),i=1..nops(gu))],guF]:

end:



#CnDnG(a,b,c,d,K): The first K terms in the sequences gcd(numer(An[i]),number(Bn[i]))
#lcm(denom(An[i]),denom(Bn[i]))
#It also returns the last 20 terms of the logs of the normailizing sequence
#as well as the last 20 terms of the implied deltas
#Try:
#CnDnG(0,0,1,0,100);
	CnDnG:=proc(a,b,c,d,K) local gu,i,lu,vu,ku,beta1,beta2:

beta1:=Lcons(c)[1]:
beta2:=1/Lcons(c)[2]:

if K<20 then
 print(K, `shpild be at least 20 `):
 RETURN(FAIL):
fi:

gu:=AnBnG(a,b,c,d,K):

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

lu:=
[
[seq(gcd(numer(gu[1][i]),numer(gu[2][i])),i=1..nops(gu[1]))],
[seq(lcm(denom(gu[1][i]),denom(gu[2][i])),i=1..nops(gu[1]))]
]:

vu:=evalf([seq((log(lu[2][i+1])-log(lu[1][i+1]))/i,i=K-20..K)]):

ku:=[ seq( (log(beta2)-vu[i])/(log(beta1)+vu[i]),i=1..nops(vu))]:

ku:=evalf(ku,10):
vu:=evalf(vu,10):

[lu,vu,ku]:

end:


#Hopeful2(): A list of generalized hopefuls
Hopefuls2:=[

[[0, 0, 1, 0]],

[[0, 0, 1/2, 0], [0, 0, 1/2, 1/2], [0, 0, 1, 1/2], [1/2, 1/2, 1/2, 0], [1/2, 1/2, 1/2, 1/2], [1/2, 1/2
, 1, 0]],
[[0, 0, 1/3, 0], [0, 0, 1/3, 1/3], [0, 0, 1/3, 2/3], [0, 0, 2/3, 0], [0, 0, 2/3, 1/3], [0, 0, 2/3, 2/3
], [0, 0, 1, 1/3], [0, 0, 1, 2/3], [0, 1/3, 1/3, 2/3], [0, 1/3, 1, 2/3], [0, 2/3, 1/3, 1/3], [0, 2/3,
1, 1/3], [1/3, 0, 1/3, 2/3], [1/3, 0, 1, 2/3], [1/3, 1/3, 1/3, 0], [1/3, 1/3, 1/3, 1/3], [1/3, 1/3, 1/
3, 2/3], [1/3, 1/3, 1, 0], [1/3, 1/3, 1, 2/3], [1/3, 2/3, 1/3, 0], [1/3, 2/3, 2/3, 0], [1/3, 2/3, 1, 0
], [2/3, 0, 1/3, 1/3], [2/3, 0, 1, 1/3], [2/3, 1/3, 1/3, 0], [2/3, 1/3, 2/3, 0], [2/3, 1/3, 1, 0], [2/
3, 2/3, 1/3, 0], [2/3, 2/3, 1/3, 1/3], [2/3, 2/3, 1/3, 2/3], [2/3, 2/3, 1, 0], [2/3, 2/3, 1, 1/3]],

[[0, 0, 1/4, 0], [0, 0, 1/4, 1/4], [0, 0, 1/4, 1/2], [0, 0, 1/4, 3/4], [0, 0, 1/2, 1/4], [0, 0, 1/2, 3
/4], [0, 0, 3/4, 1/2], [0, 0, 1, 1/4], [0, 0, 1, 3/4], [0, 1/4, 1/4, 1/2], [0, 1/2, 1/4, 0], [0, 1/2,
1/4, 1/4], [0, 1/2, 1/4, 3/4], [0, 1/2, 1/2, 1/4], [0, 1/2, 1/2, 3/4], [0, 1/2, 1, 1/4], [0, 1/2, 1, 3
/4], [0, 3/4, 1/4, 1/2], [1/4, 0, 1/4, 3/4], [1/4, 1/4, 1/4, 0], [1/4, 1/4, 1/4, 1/4], [1/4, 1/4, 1/4,
1/2], [1/4, 1/4, 1/4, 3/4], [1/4, 1/4, 1/2, 0], [1/4, 1/4, 1/2, 1/2], [1/4, 1/4, 1/2, 3/4], [1/4, 1/4,
3/4, 3/4], [1/4, 1/4, 1, 0], [1/4, 1/4, 1, 1/2], [1/4, 1/4, 1, 3/4], [1/4, 1/2, 1/4, 3/4], [1/4, 3/4,
1/4, 0], [1/4, 3/4, 1/4, 1/2], [1/4, 3/4, 1/2, 0], [1/4, 3/4, 1, 0], [1/2, 0, 1/4, 1/4], [1/2, 0, 1/4,
1/2], [1/2, 0, 1/4, 3/4], [1/2, 0, 1/2, 1/4], [1/2, 0, 1/2, 3/4], [1/2, 0, 1, 1/4], [1/2, 0, 1, 3/4],
[1/2, 1/4, 1/4, 0], [1/2, 1/2, 1/4, 0], [1/2, 1/2, 1/4, 1/4], [1/2, 1/2, 1/4, 1/2], [1/2, 1/2, 1/4, 3/
4], [1/2, 1/2, 3/4, 0], [1/2, 3/4, 1/4, 0], [3/4, 0, 1/4, 1/4], [3/4, 1/4, 1/4, 0], [3/4, 1/4, 1/4, 1/
2], [3/4, 1/4, 1/2, 0], [3/4, 1/4, 1, 0], [3/4, 1/2, 1/4, 1/4], [3/4, 3/4, 1/4, 0], [3/4, 3/4, 1/4, 1/
4], [3/4, 3/4, 1/4, 1/2], [3/4, 3/4, 1/4, 3/4], [3/4, 3/4, 1/2, 0], [3/4, 3/4, 1/2, 1/4], [3/4, 3/4, 1
/2, 1/2], [3/4, 3/4, 3/4, 1/4], [3/4, 3/4, 1, 0], [3/4, 3/4, 1, 1/4], [3/4, 3/4, 1, 1/2]],

[[0, 0, 1/5, 0], [0, 0, 1/5, 1/5], [0, 0, 1/5, 2/5], [0, 0, 1/5, 3/5], [0, 0, 1/5, 4/5], [0, 0, 2/5, 0
], [0, 0, 2/5, 1/5], [0, 0, 2/5, 2/5], [0, 0, 2/5, 3/5], [0, 0, 2/5, 4/5], [0, 0, 3/5, 0], [0, 0, 1, 1
/5], [0, 0, 1, 2/5], [0, 0, 1, 3/5], [0, 0, 1, 4/5], [0, 1/5, 1/5, 0], [0, 1/5, 1/5, 3/5], [0, 2/5, 1/
5, 0], [0, 2/5, 1/5, 1/5], [0, 3/5, 1/5, 0], [0, 3/5, 1/5, 4/5], [0, 4/5, 1/5, 0], [0, 4/5, 1/5, 2/5],
[1/5, 0, 1/5, 1/5], [1/5, 0, 1/5, 3/5], [1/5, 1/5, 1/5, 0], [1/5, 1/5, 1/5, 1/5], [1/5, 1/5, 1/5, 2/5]
, [1/5, 1/5, 1/5, 3/5], [1/5, 1/5, 1/5, 4/5], [1/5, 1/5, 2/5, 3/5], [1/5, 1/5, 2/5, 4/5], [1/5, 1/5, 1
, 3/5], [1/5, 1/5, 1, 4/5], [1/5, 2/5, 1/5, 4/5], [1/5, 2/5, 2/5, 4/5], [1/5, 2/5, 1, 4/5], [1/5, 3/5,
1/5, 2/5], [1/5, 3/5, 2/5, 2/5], [1/5, 3/5, 1, 2/5], [1/5, 4/5, 1/5, 0], [1/5, 4/5, 2/5, 0], [1/5, 4/5
, 1, 0], [2/5, 0, 1/5, 1/5], [2/5, 0, 1/5, 2/5], [2/5, 1/5, 1/5, 4/5], [2/5, 1/5, 2/5, 4/5], [2/5, 1/5
, 1, 4/5], [2/5, 2/5, 1/5, 0], [2/5, 2/5, 1/5, 1/5], [2/5, 2/5, 1/5, 2/5], [2/5, 2/5, 1/5, 3/5], [2/5,
2/5, 1/5, 4/5], [2/5, 2/5, 2/5, 1/5], [2/5, 2/5, 2/5, 3/5], [2/5, 2/5, 1, 1/5], [2/5, 2/5, 1, 3/5], [2
/5, 3/5, 1/5, 0], [2/5, 3/5, 2/5, 0], [2/5, 3/5, 1, 0], [2/5, 4/5, 1/5, 3/5], [2/5, 4/5, 2/5, 3/5], [2
/5, 4/5, 1, 3/5], [3/5, 0, 1/5, 3/5], [3/5, 0, 1/5, 4/5], [3/5, 1/5, 1/5, 2/5], [3/5, 1/5, 2/5, 2/5],
[3/5, 1/5, 1, 2/5], [3/5, 2/5, 1/5, 0], [3/5, 2/5, 2/5, 0], [3/5, 2/5, 1, 0], [3/5, 3/5, 1/5, 0], [3/5
, 3/5, 1/5, 1/5], [3/5, 3/5, 1/5, 2/5], [3/5, 3/5, 1/5, 3/5], [3/5, 3/5, 1/5, 4/5], [3/5, 3/5, 2/5, 2/
5], [3/5, 3/5, 2/5, 4/5], [3/5, 3/5, 1, 2/5], [3/5, 3/5, 1, 4/5], [3/5, 4/5, 1/5, 1/5], [3/5, 4/5, 2/5
, 1/5], [3/5, 4/5, 1, 1/5], [4/5, 0, 1/5, 2/5], [4/5, 0, 1/5, 4/5], [4/5, 1/5, 1/5, 0], [4/5, 1/5, 2/5
, 0], [4/5, 1/5, 1, 0], [4/5, 2/5, 1/5, 3/5], [4/5, 2/5, 2/5, 3/5], [4/5, 2/5, 1, 3/5], [4/5, 3/5, 1/5
, 1/5], [4/5, 3/5, 2/5, 1/5], [4/5, 3/5, 1, 1/5], [4/5, 4/5, 1/5, 0], [4/5, 4/5, 1/5, 1/5], [4/5, 4/5,
1/5, 2/5], [4/5, 4/5, 1/5, 3/5], [4/5, 4/5, 1/5, 4/5], [4/5, 4/5, 2/5, 1/5], [4/5, 4/5, 2/5, 2/5], [4/
5, 4/5, 1, 1/5], [4/5, 4/5, 1, 2/5]]]:

#Hopeful2C(r): A list of sorted generalized hopefuls
Hopefuls2C:=r -> [[[[[0, 0, 1, 0], r[1]]]], [[[[0, 0, 1, 1/2], r[1]], [[1/2, 1/2, 1, 0], 4*r[1]/(2+3*r[1])]], [[[0, 0,
1/2, 0], r[2]]], [[[0, 0, 1/2, 1/2], r[3]], [[1/2, 1/2, 1/2, 0], 8*r[3]/(4+5*r[3])]], [[[1/2, 1/2, 1/2
, 1/2], r[4]]]], [[[[0, 0, 1, 1/3], r[1]], [[0, 0, 1, 2/3], 1/2/r[1]], [[1/3, 2/3, 1, 0], -(3-6*r[1])/
(1+r[1])], [[2/3, 1/3, 1, 0], -(-6+6*r[1])/(2+r[1])]], [[[0, 0, 1/3, 0], r[2]]], [[[0, 0, 1/3, 1/3], r
[3]], [[0, 0, 1/3, 2/3], 3/4/r[3]], [[1/3, 2/3, 1/3, 0], -(27-36*r[3])/(3+2*r[3])], [[2/3, 1/3, 1/3, 0
], -(-18+18*r[3])/(2+r[3])]], [[[0, 0, 2/3, 0], r[4]]], [[[0, 0, 2/3, 1/3], r[5]], [[0, 0, 2/3, 2/3],
3/5/r[5]], [[1/3, 2/3, 2/3, 0], -(27-45*r[5])/(6+5*r[5])], [[2/3, 1/3, 2/3, 0], -(-9+9*r[5])/(2+r[5])]
], [[[0, 1/3, 1, 2/3], r[6]], [[1/3, 0, 1, 2/3], 1/2/r[6]], [[2/3, 2/3, 1, 0], -(7-14*r[6])/(2+4*r[6])
], [[2/3, 2/3, 1, 1/3], -(-2+2*r[6])/r[6]]], [[[0, 1/3, 1/3, 2/3], r[7]], [[1/3, 0, 1/3, 2/3], 3/4/r[7
]], [[2/3, 2/3, 1/3, 0], -(63-84*r[7])/(6+8*r[7])], [[2/3, 2/3, 1/3, 1/3], -(18-18*r[7])/(2-5*r[7])]],
[[[0, 2/3, 1, 1/3], r[8]], [[1/3, 1/3, 1, 0], 5*r[8]/(4+2*r[8])], [[1/3, 1/3, 1, 2/3], -3*r[8]/(-4+r[8
])], [[2/3, 0, 1, 1/3], 4*r[8]/(4+r[8])]], [[[0, 2/3, 1/3, 1/3], r[9]], [[1/3, 1/3, 1/3, 0], 15*r[9]/(
12+4*r[9])], [[1/3, 1/3, 1/3, 2/3], -3*r[9]/(-4+r[9])], [[2/3, 0, 1/3, 1/3], 12*r[9]/(12+r[9])]], [[[1
/3, 1/3, 1/3, 1/3], r[10]]], [[[2/3, 2/3, 1/3, 2/3], r[11]]]], [[[[0, 0, 1, 1/4], r[1]], [[0, 0, 1, 3/
4], 1/2/r[1]], [[1/4, 3/4, 1, 0], -(8-16*r[1])/(3+2*r[1])], [[3/4, 1/4, 1, 0], -(-8+8*r[1])/(3+r[1])]]
, [[[0, 0, 1/2, 1/4], r[2]], [[0, 0, 1/2, 3/4], 2/3/r[2]], [[1/4, 3/4, 1/2, 0], -(32-48*r[2])/(6+3*r[2
])], [[3/4, 1/4, 1/2, 0], -(-16+16*r[2])/(3+r[2])]], [[[0, 0, 1/4, 0], r[3]]], [[[0, 0, 1/4, 1/2], r[4
]], [[1/2, 1/2, 1/4, 0], 16*r[4]/(8+9*r[4])], [[1/4, 1/4, 1/4, 3/4], r[4]], [[3/4, 3/4, 1/4, 1/4], 18*
r[4]/(10+9*r[4])]], [[[0, 0, 1/4, 1/4], r[5]], [[0, 0, 1/4, 3/4], 4/5/r[5]], [[1/4, 3/4, 1/4, 0], -(
128-160*r[5])/(12+5*r[5])], [[3/4, 1/4, 1/4, 0], -(-32+32*r[5])/(3+r[5])]], [[[0, 0, 3/4, 1/2], r[6]],
[[1/2, 1/2, 3/4, 0], 16/(11+14*r[6])], [[1/4, 1/4, 3/4, 3/4], r[6]], [[3/4, 3/4, 3/4, 1/4], -(-6+3*r[6
])/(2+4*r[6])]], [[[0, 1/2, 1, 1/4], r[7]], [[1/2, 0, 1, 1/4], 5*r[7]/(5+r[7])], [[1/4, 1/4, 1, 0], 6*
r[7]/(5+2*r[7])], [[1/4, 1/4, 1, 1/2], -4*r[7]/(-5+r[7])]], [[[0, 1/2, 1, 3/4], r[8]], [[1/2, 0, 1, 3/
4], 1/2/r[8]], [[3/4, 3/4, 1, 0], -(10-20*r[8])/(3+6*r[8])], [[3/4, 3/4, 1, 1/2], -(8-8*r[8])/(1-5*r[8
])]], [[[0, 1/2, 1/2, 1/4], r[9]], [[1/2, 0, 1/2, 1/4], 10*r[9]/(10+r[9])], [[1/4, 1/4, 1/2, 0], 12*r[
9]/(10+3*r[9])], [[1/4, 1/4, 1/2, 1/2], -4*r[9]/(-5+r[9])]], [[[0, 1/2, 1/2, 3/4], r[10]], [[1/2, 0, 1
/2, 3/4], 2/3/r[10]], [[3/4, 3/4, 1/2, 0], -(40-60*r[10])/(6+9*r[10])], [[3/4, 3/4, 1/2, 1/2], -(16-16
*r[10])/(3-7*r[10])]], [[[0, 1/2, 1/4, 0], r[11]], [[1/2, 0, 1/4, 1/2], -2*r[11]/(-3+r[11])]], [[[0, 1
/2, 1/4, 1/4], r[12]], [[1/2, 0, 1/4, 1/4], 20*r[12]/(20+r[12])], [[1/4, 1/4, 1/4, 0], 24*r[12]/(20+5*
r[12])], [[1/4, 1/4, 1/4, 1/2], -4*r[12]/(-5+r[12])]], [[[0, 1/2, 1/4, 3/4], r[13]], [[1/2, 0, 1/4, 3/
4], 4/5/r[13]], [[3/4, 3/4, 1/4, 0], -(160-200*r[13])/(12+15*r[13])], [[3/4, 3/4, 1/4, 1/2], -(32-32*r
[13])/(7-11*r[13])]], [[[0, 1/4, 1/4, 1/2], r[14]], [[1/2, 3/4, 1/4, 0], -(72-90*r[14])/(6+5*r[14])],
[[1/4, 0, 1/4, 3/4], 4/5/r[14]], [[3/4, 1/2, 1/4, 1/4], -(-16+16*r[14])/(-1+3*r[14])]], [[[0, 3/4, 1/4
, 1/2], r[15]], [[1/2, 1/4, 1/4, 0], 14*r[15]/(10+5*r[15])], [[1/4, 1/2, 1/4, 3/4], -4*r[15]/(-5+r[15]
)], [[3/4, 0, 1/4, 1/4], 12*r[15]/(10+3*r[15])]], [[[1/2, 1/2, 1/4, 1/2], r[16]]], [[[1/2, 1/2, 1/4, 1
/4], r[17]], [[1/2, 1/2, 1/4, 3/4], -20*r[17]/(-28+9*r[17])], [[1/4, 3/4, 1/4, 1/2], -24*r[17]/(-28+5*
r[17])], [[3/4, 1/4, 1/4, 1/2], -6*r[17]/(-7+r[17])]], [[[1/4, 1/4, 1, 3/4], r[18]], [[3/4, 3/4, 1, 1/
4], 9*r[18]/(5+6*r[18])]], [[[1/4, 1/4, 1/2, 3/4], r[19]], [[3/4, 3/4, 1/2, 1/4], 9*r[19]/(5+5*r[19])]
], [[[1/4, 1/4, 1/4, 1/4], r[20]]], [[[3/4, 3/4, 1/4, 3/4], r[21]]]], [[[[0, 0, 1, 1/5], r[1]], [[0, 0
, 1, 4/5], 1/2/r[1]], [[1/5, 4/5, 1, 0], -(5-10*r[1])/(2+r[1])], [[4/5, 1/5, 1, 0], -(-10+10*r[1])/(4+
r[1])]], [[[0, 0, 1, 2/5], r[2]], [[0, 0, 1, 3/5], 1/2/r[2]], [[2/5, 3/5, 1, 0], -(10-20*r[2])/(3+4*r[
2])], [[3/5, 2/5, 1, 0], -(-10+10*r[2])/(3+2*r[2])]], [[[0, 0, 1/5, 0], r[3]]], [[[0, 0, 1/5, 1/5], r[
4]], [[0, 0, 1/5, 4/5], 5/6/r[4]], [[1/5, 4/5, 1/5, 0], -(125-150*r[4])/(10+3*r[4])], [[4/5, 1/5, 1/5,
0], -(-50+50*r[4])/(4+r[4])]], [[[0, 0, 1/5, 2/5], r[5]], [[0, 0, 1/5, 3/5], 5/6/r[5]], [[2/5, 3/5, 1/
5, 0], -(250-300*r[5])/(15+12*r[5])], [[3/5, 2/5, 1/5, 0], -(-50+50*r[5])/(3+2*r[5])]], [[[0, 0, 2/5,
0], r[6]]], [[[0, 0, 2/5, 1/5], r[7]], [[0, 0, 2/5, 4/5], 5/7/r[7]], [[1/5, 4/5, 2/5, 0], -(125-175*r[
7])/(20+7*r[7])], [[4/5, 1/5, 2/5, 0], -(-25+25*r[7])/(4+r[7])]], [[[0, 0, 2/5, 2/5], r[8]], [[0, 0, 2
/5, 3/5], 5/7/r[8]], [[2/5, 3/5, 2/5, 0], -(125-175*r[8])/(15+14*r[8])], [[3/5, 2/5, 2/5, 0], -(-25+25
*r[8])/(3+2*r[8])]], [[[0, 0, 3/5, 0], r[9]]], [[[0, 1/5, 1/5, 0], r[10]], [[1/5, 0, 1/5, 1/5], -5*r[
10]/(-6+r[10])]], [[[0, 1/5, 1/5, 3/5], r[11]], [[1/5, 0, 1/5, 3/5], 5/6/r[11]], [[3/5, 3/5, 1/5, 0],
-(275-330*r[11])/(15+18*r[11])], [[3/5, 3/5, 1/5, 1/5], -(-50+50*r[11])/(-3+8*r[11])]], [[[0, 2/5, 1/5
, 0], r[12]], [[2/5, 0, 1/5, 2/5], -5*r[12]/(-7+2*r[12])]], [[[0, 2/5, 1/5, 1/5], r[13]], [[1/5, 1/5,
1/5, 0], 35*r[13]/(30+6*r[13])], [[1/5, 1/5, 1/5, 2/5], -5*r[13]/(-6+r[13])], [[2/5, 0, 1/5, 1/5], 30*
r[13]/(30+r[13])]], [[[0, 3/5, 1/5, 0], r[14]], [[3/5, 0, 1/5, 3/5], -5*r[14]/(-8+3*r[14])]], [[[0, 3/
5, 1/5, 4/5], r[15]], [[3/5, 0, 1/5, 4/5], 5/6/r[15]], [[4/5, 4/5, 1/5, 0], -(325-390*r[15])/(20+24*r[
15])], [[4/5, 4/5, 1/5, 3/5], -(50-50*r[15])/(14-19*r[15])]], [[[0, 4/5, 1/5, 0], r[16]], [[4/5, 0, 1/
5, 4/5], -5*r[16]/(-9+4*r[16])]], [[[0, 4/5, 1/5, 2/5], r[17]], [[2/5, 2/5, 1/5, 0], 45*r[17]/(35+12*r
[17])], [[2/5, 2/5, 1/5, 4/5], -5*r[17]/(-7+2*r[17])], [[4/5, 0, 1/5, 2/5], 35*r[17]/(35+2*r[17])]], [
[[1/5, 1/5, 1, 3/5], r[18]], [[1/5, 1/5, 1, 4/5], 1/2/r[18]], [[3/5, 4/5, 1, 1/5], -(11-22*r[18])/(5+2
*r[18])], [[4/5, 3/5, 1, 1/5], -(-11+11*r[18])/(2+4*r[18])]], [[[1/5, 1/5, 1/5, 1/5], r[19]]], [[[1/5,
1/5, 1/5, 3/5], r[20]], [[1/5, 1/5, 1/5, 4/5], 5/6/r[20]], [[3/5, 4/5, 1/5, 1/5], -(275-330*r[20])/(45
-18*r[20])], [[4/5, 3/5, 1/5, 1/5], -(55-55*r[20])/(2-8*r[20])]], [[[1/5, 1/5, 2/5, 3/5], r[21]], [[1/
5, 1/5, 2/5, 4/5], 5/7/r[21]], [[3/5, 4/5, 2/5, 1/5], -(275-385*r[21])/(65-7*r[21])], [[4/5, 3/5, 2/5,
1/5], -(-55+55*r[21])/(1+11*r[21])]], [[[1/5, 2/5, 1, 4/5], r[22]], [[2/5, 1/5, 1, 4/5], 1/2/r[22]], [
[4/5, 4/5, 1, 1/5], -(12-24*r[22])/(5+4*r[22])], [[4/5, 4/5, 1, 2/5], -1/6*(-11+11*r[22])/r[22]]], [[[
1/5, 2/5, 1/5, 4/5], r[23]], [[2/5, 1/5, 1/5, 4/5], 5/6/r[23]], [[4/5, 4/5, 1/5, 1/5], -(-100+120*r[23
])/(-15+4*r[23])], [[4/5, 4/5, 1/5, 2/5], -(-55+55*r[23])/(-8+14*r[23])]], [[[1/5, 2/5, 2/5, 4/5], r[
24]], [[2/5, 1/5, 2/5, 4/5], 5/7/r[24]], [[4/5, 4/5, 2/5, 1/5], -(300-420*r[24])/(65+7*r[24])], [[4/5,
4/5, 2/5, 2/5], -(-55+55*r[24])/(-6+18*r[24])]], [[[1/5, 3/5, 1, 2/5], r[25]], [[2/5, 2/5, 1, 1/5], 8*
r[25]/(7+2*r[25])], [[2/5, 2/5, 1, 3/5], -6*r[25]/(-7+r[25])], [[3/5, 1/5, 1, 2/5], 7*r[25]/(7+r[25])]
], [[[1/5, 3/5, 1/5, 2/5], r[26]], [[2/5, 2/5, 1/5, 1/5], 40*r[26]/(35+6*r[26])], [[2/5, 2/5, 1/5, 3/5
], -6*r[26]/(-7+r[26])], [[3/5, 1/5, 1/5, 2/5], 35*r[26]/(35+r[26])]], [[[1/5, 3/5, 2/5, 2/5], r[27]],
[[2/5, 2/5, 2/5, 1/5], 40*r[27]/(35+7*r[27])], [[2/5, 2/5, 2/5, 3/5], -6*r[27]/(-7+r[27])], [[3/5, 1/5
, 2/5, 2/5], 35*r[27]/(35+2*r[27])]], [[[2/5, 2/5, 1/5, 2/5], r[28]]], [[[2/5, 4/5, 1, 3/5], r[29]], [
[3/5, 3/5, 1, 2/5], 9*r[29]/(8+2*r[29])], [[3/5, 3/5, 1, 4/5], -7*r[29]/(-8+r[29])], [[4/5, 2/5, 1, 3/
5], 8*r[29]/(8+r[29])]], [[[2/5, 4/5, 1/5, 3/5], r[30]], [[3/5, 3/5, 1/5, 2/5], 45*r[30]/(40+6*r[30])]
, [[3/5, 3/5, 1/5, 4/5], -7*r[30]/(-8+r[30])], [[4/5, 2/5, 1/5, 3/5], 40*r[30]/(40+r[30])]], [[[2/5, 4
/5, 2/5, 3/5], r[31]], [[3/5, 3/5, 2/5, 2/5], 45*r[31]/(40+7*r[31])], [[3/5, 3/5, 2/5, 4/5], -7*r[31]/
(-8+r[31])], [[4/5, 2/5, 2/5, 3/5], 20*r[31]/(20+r[31])]], [[[3/5, 3/5, 1/5, 3/5], r[32]]], [[[4/5, 4/
5, 1/5, 4/5], r[33]]]]]:



#SortCG(L,N,r): Given a list of quadruples, L, a positive integer N, and a symbol r
#divides them into classes such that GLCg(op(L[i]))[1] are related to
#each other (conjecturally) by a fractional-linear transfomation with integers  less than N in absolute value
#The output is a list of lists such that the first entry is the main quadruple, followed by the symbol r[i],
#its decimal approximation, followed by the tuples, and the expression in terms of r[i]
#Try:
#SortCG(Hopefuls2()[2],10000,r);
SortCG:=proc(L, N, r)
    local i,S1,Rf,lu1,lu2,co,T,ka:
    S1:=convert(L,set):

    co:=0:
    while S1<>{} do
    co:=co+1:
     lu1:=S1[1]:
     Rf:=GLCg(op(lu1)):
     T[co]:=[]:

     for lu2 in S1   do
       ka:=MyIDs(GLCg(op(lu2)),r[co],Rf,N):
      if ka<>FAIL then
       T[co]:=[op(T[co]),[lu2,ka]]:
      S1:=S1 minus {lu2}:
     fi:
    od:

    od:

    [seq(T[i],i=1..co)]:

end:


#Hopefuls2v: A description of generalized hopefuls. Try:
#Hopefuls2v():

Hopefuls2v:=proc() local r,gu,C,lu,i,j,con1,con2,S1,S2:

print(`A list of Hopeful quartets for Int(x^a*(1-x)^b/(1+c*x)^(d+1),x=0..1)/Int(x^a*(1-x)^b/(1+c*x)^(d),x=0..1)`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Definition: For rational a,b,c,d, let `):
print(`C(a,b,c,d)=Int(x^a*(1-x)^b/(1+c*x)^(d+1),x=0..1)/Int(x^a*(1-x)^b/(1+c*x)^d,x=0..1)`):

gu:=Hopefuls2C(r):

S1:=[]:
S2:=[]:

for i from 1 to nops(gu) do
print(``):
print(`There are `, nops(gu[i]), `inequivalent classes with denominator `, i):

print(`Here there are`):

for j from 1 to nops(gu[i]) do

 lu:=gu[i][j][1][1]:


con1:=GLCg(op(lu)):

con2:=identify(con1):


if con1<>con2 then
 print(C(op(lu))=con2):
S1:=[op(S1),[C(op(lu)),con2]]:
else
 print(C(op(lu))):
S2:=[op(S2),C(op(lu))]:
fi:
od:
od:

print(``):

print(`To sum up, here are the identified 4-tuples`):
print(``):
lprint(S1):
print(``):
print(` here are the unidentified 4-tuples, that are candidates to first-ever irrationality proof`):
print(``):
lprint(S2):
print(``):
end: