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

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



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

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

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


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

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

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

with(combinat):



ezraSt:=proc()

if args=NULL then
 print(` The  STORY procedures are: PaperA, PaperQ, PaperAT, PaperATS `):
 print(``):

else
ezra(args):
fi:

end:


ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AZd1, Bas, Check37, CutO, CutOa, CutOb, DeltSp, Frat, GuessK, GuessPK, Hope, TripSeqSp, TupSeq, ZugSeq `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: APx , AppSeq, APxF, APxFS, AZd , Delta1S, Delta2S, deltE, deltP, deltT , FindGoodP,FindGoodPAT, FindGoodPg `):
 print(` `):



elif nops([args])=1 and op(1,[args])=AppSeq then
print(`AppSeq(P,x,N): given a polynomial P in x of degree 1 or 2, output c= APx(P,x,0) (the irrational number we wish to approximate), `):
print(`and the sequence of diophantine approximations inspired by APx(P,x,n); Try:`):
print(`AppSeq(1+x,x,30);`):

elif nops([args])=1 and op(1,[args])=APx then
print(`APx(P,x,n1): Given a polynomial P of x , outputs int(x^n1*(1-x)^n1/P(x)^(n1+1),x=0..1); Done directly. Try:`):
print(`[seq(APx(1+x,x,n1),n1=0..10)];`):



elif nops([args])=1 and op(1,[args])=APxF then
print(`APxF(P,x,n1): Given a polynomial P of x , outputs int(x^n1*(1-x)^n1/P(x)^(n1+1),x=0..1); Done via the recurrence obtained `):
print(`From the Alkvist-Zeilberger algorithm. Try:`):
print(`[seq(APxF(1+x,x,n1),n1=0..10)];`):


elif nops([args])=1 and op(1,[args])=APxFS then 
print(`APxFS(P,x,n1,X): Like APxF(P,x,n1) but in terms of the initial conditions X[0], ..., X[ORD-1] where ORD is the order of the recurrence`):
print(`Try:`):
print(`APxFS(1+x+x^2+x^3,x,10,X);`):

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 nops([args])=1 and op(1,[args])=AZd1 then
print(`AZd1(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 ope(n,N^(-1)) such that`):
print(`The f(n):=Int( A(y,n),y); (assuming a contour integral or one that vanishes at the end points`):
print(`satisfies: f(n)=ope(n,N^(-1))f(n). Try:`):
print(`AZd1(x^n*(1-x)^n,x,n,N);`):


elif nops([args])=1 and op(1,[args])=Bas then
print(`Bas(P,x): inputs a cubic polynomial in x with integer coefficients and outputs two constants {theta1,theta2} such that`):
print(`APx(F,x,n1) is always a linear combination of {1,theta1,theta2} with rational coefficients. Try:`):
print(`Bas(1+x+x^2+x^3,x);`):


elif nops([args])=1 and op(1,[args])=Check37 then
print(`Check37(L): checks that the conjecture that for all integers a that are 1,3,7 mod 8 up to L`):
print(`the divisibility constant of x^2+a is 2*sqrt(a). It outputs the exceptios. It the output`):
print(`is the empty set it is good news. Try:`):
print(`Check37(10);`):

elif nops([args])=1 and op(1,[args])=CutO then
print(`CutO(a1): Inputs a and outputs the largest b such that Delta1S(a,b) is positive. Try:`):
print(`CutO(5);`):

elif nops([args])=1 and op(1,[args])=CutOa then
print(`CutOa(b): the smallest a such that Delta1S(a,b) is positive. Try:`):
print(`CutOa(b);`):

elif nops([args])=1 and op(1,[args])=CutOb then
print(`CutOb(a): the smallest a such that Delta1S(a,b) is positive. Try:`):
print(`CutOb(a);`):


elif nops([args])=1 and op(1,[args])=Delta1S then
print(`Delta1S(a,b):The explicit expression, in terms of a and b,  for the delta in the approximation of`):
print(`log((a+b)/a) implied by the integral int(x^n*(1-x)^n)/(a+b*x)^(n+1),x=0..1). Try:`):
print(`Delta1S(a,b); `):


elif nops([args])=1 and op(1,[args])=Delta2S then
print(`Delta2S(a,b,c,K):The explicit expression, in terms of a and b,  for the delta in the approximation of`):
print(` int(1/(a+b*x+c*x^2),x=0..1), implied by the integral int(x^n*(1-x)^n)/(a+b*x+c*x^2)^(n+1),x=0..1). `):
print(`assuming that the divisibility constant is K. This is SYMBOLIC. To implement it you need to know K.`):
print(`Delta2S(a,b,c,K);`):

elif nops([args])=1 and op(1,[args])=DeltSp then
print(`DeltSp(): The exponent such that APxF(1+x+x^2+x^3,x) gives a simultaneous approximation to 1, log(2), Pi. Try:`):
print(`DeltSp();`):

elif nops([args])=1 and op(1,[args])=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 nops([args])=1 and op(1,[args])=deltP  then
print(`The sequence of deltas for AppSeq from N-20 to N. Try:`):
print(`deltP(1+x,x,50);`):


elif nops([args])=1 and op(1,[args])=deltT then
print(`deltT(P,x): inputs a polynomial P of degree 1 or 2, finds the theoretical delta for   the constant int(1/P,x=0..1) implied by`):
print(`the corresponding sequence int(x^n*(1-x)^n/P^(n+1),x=0..1); It returns, c, the conjectured K, and the delta, and its floating-point.`):
print(`So the output is`):
print(`[c,K,delta,evalf(delta), SmallestEmpricalDelta(from 180 to 200)]: Try:`):
print(`deltT(1+x,x);`):

elif nops([args])=1 and op(1,[args])=Frat then
print(`Frat(L): given a linear combination of constants and a rational number, exctracts the rational number. Try:`):
print(`Frat(1/11+log(5)+log(11)+Pi);`):



elif nops([args])=1 and op(1,[args])=FindGoodP then
print(`FindGoodP(L,x): ouputs the list of [a+b*x+c*x^2,constant, K, ExactValueOfDelta, ImpliedIrrationalityMeasure] such that the sequence `):
print(` int( x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant `):
print(`int( 1/(a+b*x+c*x^2),x=0..1) for all 1<=a,b,c <= L. We only list irreducible quardratics, since otherwise it is covered by PaperA(). Try: `):
print(`FindGoodP(2,x);`):


elif nops([args])=1 and op(1,[args])=FindGoodPAT then
print(`FindGoodPAT(L,x): ouputs the list of [a+c*x^2,constant, K, ExactValueOfDelta, ImpliedIrrationalityMeasure] such that the sequence `):
print(` int( x^n*(1-x)^n/(a+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant `):
print(`int( 1/(a+c*x^2),x=0..1) (expressible in terms of arctan, hence the name of the procedure)`):
print(`for all 1<=a,b,c <= L. Try: `):
print(`FindGoodPAT(3,x);`):


elif nops([args])=1 and op(1,[args])=FindGoodPg then
print(`FindGoodPg(L): ouputs the list of [a+b*x+c*x^2,K,delta] such that the sequence `):
print(` int( x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant `):
print(` int( 1/(a+b*x+c*x^2),x=0..1) for all 1<=a<=L, c<>0, and -L<=b,c <= L. Try: `):
print(`FindGoodPg(2);`):

elif nops([args])=1 and op(1,[args])=GuessPK then
print(`GuessPK(P,x,N): Given a polynomial P in x of degree 1 or 2 and a positive integer N (at least 40) guesses the`):
print(`constant such that the coefficients of int(x^n*(1-x)^n/P^(n+1),x=0..1) when multiplied by lcm(1...n)*K^n integers. Try:`):
print(`GuessPK(1+x,x,40);`):
print(`GuessPK(1+x+x^2,x,40);`):

elif nops([args])=1 and op(1,[args])=GuessK then
print(`GuessK(Lis): Inputs a list of rational numbers, Lis, and outputs a positive integer K, and a small integer L such that`):
print(`Lis[n]*lcm(1...n)*K^n*L are all integers. Try: `):
print(`GuessK([seq(5^i/(dn(i)*7), i=1..30)]); `):

elif nops([args])=1 and op(1,[args])=Hope then
print(`Hope(P,x,N): inputs a polynomial P of x with integer coefficients of degree 1 or 2 in x and a positive integer`):
print(`outputs the empircal delta for the N-th approximation. If it is positive there is hope that `):
print(`we have an irratioality proof for the given constant. `):
print(`The output is the constant and the empircal delta. Try:`):
print(`Hope(1+x,x,100);`):
print(`Hope(13+7*x+x^2,x,100);`):



elif nops([args])=1 and op(1,[args])=PaperA then
print(`PaperA(): A computer-generated paper about the irrationality of log(1+b/a) for integer b< (2*sqrt(a*e)+1)/e`):
print(`and an explicit irratinality measure. It is a computer-generated version of Theorem 1 in  `):
print(`K. Alladi and M.L.Robinson, Legendre polynomials and irrationality, Crelle J.  318 (1980), 134-155.`):
print(`but is more stream-lined, direct and does NOT mention  Legedre polynomials. Type:`):
print(`PaperA(); `):


elif nops([args])=1 and op(1,[args])=PaperAT then
print(`PaperAT(L): Inputs a positive integer L, outputs an article about the Diophantine approximations of Int(1/(a+c*x^2),x=0..1) for 1<=a,c<=L, and`):
print(` that yield irrationality proofs and their corresponding measures. Try: `):
print(`PaperAT(3);`):

elif nops([args])=1 and op(1,[args])=PaperATS then
print(`PaperATS(): A computer-generated paper about the irrationality of `):
print(`arctan(sqrt(a))/ sqrt(a), alias, int(1/(x^2+a),x=0..1) for ALL positive integers a that are 3 mod 8 or 7 mod 8.`):
print(`It gives a symbolic expression in a for the irrationality measure, complete with proof`):
print(`(modulo a divisibility lemma left to the reader). Try:`):
print(`PaperATS();`):

elif nops([args])=1 and op(1,[args])=PaperQ then
print(`PaperQ(L): Inputs a positive integer L, outputs an article about the Diophantine approximations of Int(1/(a+b*x+c*x^2),x=0..1) for 1<=a,b,c<=L, and`):
print(` that yield irrationality proofs and their corresponding measures. Try: `):
print(`PaperQ(3);`):


elif nops([args])=1 and op(1,[args])=TripSeqSp then
print(`TripSeqSp(N): the sequence of coefficients of 1, log(2), and Pi for APx(1+x+x^2+x^3,x,n) for n from 1 to N. Try:`):
print(`TripSeqSp(30);`):


elif nops([args])=1 and op(1,[args])=TupSeq then
print(`TupSeq(P,x,N): given a polynomial P in x  output the sequence of coefficients of X[0], .., X[Order-1] in`):
print(`APxFS(P,x,n) from n=1 to N. Try:`):
print(`TupSeq(1+x+x^2+x^3,x,30);`):

elif nops([args])=1 and op(1,[args])=ZugSeq then
print(`ZugSeq(P,x,N): given a polynomial P in x of degree 1 or 2, output sAPx(P,x,0) (the irrational number we wish to approximate), let't call it c`):
print(`and the sequence of pairs [An,Bn] from n=1 to N such that APx(P,x,n)=An+Bn*c. `):
print(`The output is the pair [c,SequenceOfPairs]. Try: `):
print(`ZugSeq(1+x,x,30);`):


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

Digits:=7000:

#########From EKHAD

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

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

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


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

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

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

0:
end:

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


dn:=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:

 
#APx(P,x,n1): Given a polynomial P of x , outputs int(x^n1*(1-x)^n1/P(x)^(n1+1),x=0..1); Try:
#APx(1+x,x,5);
APx:=proc(P,x,n1)
int(x^n1*(1-x)^n1/P^(n1+1),x=0..1):
end:

 
#APxF(P,x,n1): A fast version of APx(P,x,n1) using a recurrence obtatained with the Almkvist-Zeilberger algorithm. Try:
#[seq(APxF(1+x,x,n1),n1=0..20)]:
APxF:=proc(P,x,n1) local ope,N,ORD1,n,i1:
option remember:
ope:=AZd1(x^n*(1-x)^n/P^(n+1),x,n,N):
ORD1:=-ldegree(ope,N):

if n1<=ORD1-1 then
 APx(P,x,n1):
else
add(subs(n=n1,coeff(ope,N,-i1))*APxF(P,x,n1-i1),i1=1..ORD1):
fi:

end:

#Frat(L): given a linear combination of constants and a rational number, exctracts the rational number. Try:
#Frat(1/11+log(5)+log(11)+Pi);
Frat:=proc(L) local i,L1:

L1:=expand(L):

if ( type(L1,fraction) or type(L1, integer) ) then
 RETURN(L1):
fi:

if type(L1,`*`) then
 RETURN(0):
fi:


if type(L1,function) then
 RETURN(0):
fi:

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

for i from 1 to nops(L1) do
 if (type(op(i,L1),fraction)  or type(op(i,L1),integer) ) then
  RETURN(op(i,L1)):
 fi:
od:

0:
end:


#ZugSeq(P,x,N): given a polynomial P in x of degree 1 or 2, output sAPx(P,x,0) (the irrational number we wish to approximate), let't call it c
#and the sequence of pairs [An,Bn] from n=1 to N such that APx(P,x,n)=An+Bn*c. Try:
#ZugSeq(1+x,x,30);
ZugSeq:=proc(P,x,N) local c,n1,lu,mu1,mu2,gu:

if not (degree(P,x)=1 or degree(P,x)=2) then
 print(P, `should be a polynomial in`, x, `of degree 1 or 2 `):
 RETURN(FAIL):
fi:
c:=APx(P,x,0):
gu:=[]:

for n1 from 1 to N do
 lu:=APxF(P,x,n1):
 mu1:=Frat(lu):
 mu2:=simplify((lu-mu1)/c):

 if not (type(mu2,fraction) or type(mu2,integer)) then
  RETURN(FAIL):
 fi:

 if simplify(lu-mu1-mu2*c)<>0 then
    RETURN(FAIL):
fi:

 gu:=[op(gu),[mu1,mu2]]:
od:

[c,gu]:

end:



#AppSeq(P,x,N): given a polynomial P in x of degree 1 or 2, output c= APx(P,x,0) (the irrational number we wish to approximate), 
#and the sequence of diophantine approximations inspired by APx(P,x,n); Try:
#AppSeq,30);
AppSeq:=proc(P,x,N) local gu,c, i:

gu:=ZugSeq(P,x,N):

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

c:=gu[1]:

gu:=gu[2]:

gu:=[seq(-gu[i][1]/gu[i][2],i=1..N)]:
[c,gu]:
end:


#The sequence of deltas for AppSeq from N-20 to N. Try:
#deltP(1+x,x,50);
deltP:=proc(P,x,N) local gu,c,i:


gu:=AppSeq(P,x,N):

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

if N<=20 then
 print(N, `must be at least 20 `):
fi:

c:=gu[1]:

gu:=gu[2]:

[seq(deltE(gu[i],c),i=N-20..N)]:

end:

#Hope(P,x,N): inputs a polynomial P of x with integer coefficients of degree 1 or 2 in x and a positive integer
#outputs the empircal delta for the N-th approximation. If it is positive there is hope that 
#we have an irratioality proof for the given constant. 
#The output is the constant and the empircal delta. Try:
#Hope(1+x,x,100);
Hope:=proc(P,x,N) local gu,c:

gu:=AppSeq(P,x,N):

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

c:=gu[1]:
gu:=gu[2]:

[c,deltE(gu[N],c)]:

end:

#GuessK1(Lis): Inputs a list of rational numbers, Lis, and outputs a positive integer K, and a small integer L such that
#Lis[n]*lcm(1...n)*K^n*L are all integers. Try
#GuessK1([seq(5^i/(dn(i)*7), i=1..30)]);
GuessK1:=proc(Lis) local N,n1,gu,K,i,mu,vu:

N:=nops(Lis):
if N<=21 then
 print(`The list has to be at least of length 21`):
fi:

gu:=ifactors(denom(Lis[N]*dn(N)))[2]:


mu:=1:

for i from 1 to nops(gu) do
 vu:=round(gu[i][2]/N):

  if abs(vu-gu[i][2]/N)<1/2 then
   mu:=mu*gu[i][1]^vu:
  fi:
od:
 

K:=mu:


gu:=lcm(seq(denom(Lis[n1]*dn(n1)*K^n1),n1=20..N)):


if {seq(denom(gu*Lis[n1]*dn(n1)*K^n1),n1=20..N)}<>{1} then
 RETURN(FAIL):
fi:

[K, gu]:


end:



#GuessK(Lis): Inputs a list of rational numbers, Lis, and outputs a positive integer K, and a small integer L such that
#Lis[n]*lcm(1...n)*K^n*L are all integers. Try
#GuessK([seq(5^i/(dn(i)*7), i=1..30)]);
GuessK:=proc(Lis) local N,K,gu,i,lu,ru:
K:=GuessK1(Lis):

if K[2]<=1000 then

if max(seq(denom(K[1]^i*Lis[i]*dn(i)*K[2]),i=1..nops(Lis)))>1000 then
  RETURN(FAIL):
 else
 RETURN(K):
fi:

fi:

N:=nops(Lis):

gu:=ifactors(denom(Lis[N]*dn(N)*K[1]^N))[2]:



lu:=1:

for i from 1 to nops(gu) do
 if abs(gu[i][2]/N-1/2)<1/5 then
  lu:=lu*gu[i][1]:
 fi:
od:


ru:=max(seq(denom(Lis[i]*dn(i)*K[1]^i*lu^(trunc(i/2))),i=1..nops(Lis))):



[K[1]*sqrt(lu),ru]:
end:








#deltT(P,x): inputs a polynomial P of degree 1 or 2, finds the theoretical delta for   the constant int(1/P,x=0..1) implied by
#the corresponding sequence int(x^n*(1-x)^n/P^(n+1),x=0..1); It returns, c, the conjectured K, and the delta, and its floating-point.
#So the output is
#[c,K,delta,evalf(delta), ), SmallestEmpricalDelta(from 180 to 200) ]: Try:
#deltT(1+x,x);
deltT:=proc(P,x) local gu,K1,K2,K,ope,c,N,i,n,lu:

gu:=ZugSeq(P,x,100):
c:=gu[1]:
gu:=gu[2]:



K1:=GuessK([seq(gu[i][1]*dn(i),i=1..nops(gu))]):

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

K2:=GuessK([seq(gu[i][2]*dn(i),i=1..nops(gu))]):

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




if K1[2]>1000 or K2[2]>1000 then
 RETURN(FAIL):
fi:

if K1[1]<>K2[1] then
 print(`The K-s are not the same`):
fi:

if type(K1[1],integer) and type(K2[1],integer) then
 K:=lcm(K1[1],K2[1]):
elif type(K1[1]^2,integer) and type(K2[1]^2,integer) then
 K:=sqrt(lcm(K1[1]^2,K2[1]^2)):
else
 RETURN(FAIL):
fi:




ope:=AZd(x^n*(1-x)^n/P^(n+1),x,n,N)[1]:




if degree(ope,N)<>2 then
 print(`The operator`, ope, `is not second-order `):
 RETURN(FAIL):
fi:

ope:=lcoeff(ope,n):



gu:=[solve(ope,N)]:

gu:=[abs(gu[1]), abs(gu[2])]:

if evalf(gu[1])>evalf(gu[2]) then
 gu:=[gu[2],gu[1]]:
fi:


lu:=-(log(gu[1])+log(K)+1)/(log(gu[2])+log(K)+1):

[c,K,lu,evalf(lu,20), min(op(deltP(P,x,100))) ]:

end:





#Delta1S(a,b):The explicit expression, in terms of a and b,  for the delta in the approximation of
#log((a+b)/a) implied by the integral int(x^n*(1-x)^n)/(a+b*x)^(n+1),x=0..1). Try:
#Delta1S(a,b);
Delta1S:=proc(a,b) local ope,gu,gu1,lu,x,n,N:

ope:=AZd((x^n*(1-x)^n)/(a+b*x)^(n+1),x,n,N)[1];

ope:=lcoeff(ope,n):

gu:=[solve(ope,N)]:


gu1:=subs({a=2,b=3},gu):

if evalf(abs(gu1[1]))>evalf(abs(gu1[2])) then
 gu:=[gu[2],gu[1]]:
fi:

lu:=-(log(gu[1])+2*log(b)+1)/(log(gu[2])+2*log(b)+1):

lu:
end:




#Delta2S(a,b,c,K):The explicit expression, in terms of a and b,  for the delta in the approximation of
# int(1/(a+b*x+c*x^2),x=0..1), implied by the integral int(x^n*(1-x)^n)/(a+b*x+c*x^2)^(n+1),x=0..1). 
#assuming that the divisibility constant is K. This is SYMBOLIC. To implement it you need to know K.
#Delta2S(a,b,c,K);
Delta2S:=proc(a,b,c,K) local ope,gu,gu1,lu,x,n,N:

ope:=AZd((x^n*(1-x)^n)/(a+b*x+c*x^2)^(n+1),x,n,N)[1];

ope:=lcoeff(ope,n):

gu:=[solve(ope,N)]:

gu:=[abs(gu[1]),abs(gu[2])]:

gu1:=subs({a=2,b=3,c=4},gu):

if evalf(gu1[1])>evalf(gu1[2]) then
 gu:=[gu[2],gu[1]]:
fi:

lu:=-(log(gu[1])+log(K)+1)/(log(gu[2])+log(K)+1):

lu:
end:


#CutO(a1): Inputs a and outputs the largest b such that Delta1S(a,b) is positive. Try:
#CutO(5);
CutO:=proc(a1) local lu,a,b,b1:
 lu:=Delta1S(a,b):

  for b1 from 1 while evalf(subs({a=a1,b=b1},lu))>0 do od:
b1-1:
end:

#CutOa(b): the smallest a such that Delta1S(a,b) is positive. Try:
#CutOa(b);
CutOa:=proc(b): ((b-exp(-1))/2)^2:end:

#CutOb(a): the largest b such that Delta1S(a,b) is positive. Try:
#CutOb(a);
CutOb:=proc(a) : 
(2*sqrt(a*exp(1))+1)/exp(1):
end:


#PaperA(): A computer-generated paper about the irrationality of log(1+b/a) for integer b< (2*sqrt(a*e)+1)/e
#and an explicit irratinality measure. It is a computer-generated version of Theorem 1 in  
#K. Alladi and M.L.Robinson, Legendre polynomials and irrationality, Crelle J.  318 (1980), 134-155.
#but is more stream-lined, direct and does NOT mention  Legedre polynomials. Type:
#PaperA();
PaperA:=proc() local E,n,a,b,de,ope,N,i,x,gu,gu1:
print(`Irrationality Proofs and Measures for log(1+b/a) for positive integers a and b such that a>(b-1/e)^2/4 `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

de:=Delta1S(a,b):

print(`Theorem: Let a and b be positive integers such that a>(b-1/e)^2/4, then log(1+b/a) is an irrational number`):
print(`and an irrationality measure is`, simplify(1+1/de) ):
print(``):
print(`In Maple format the irrationality measure is`):
print(``):
lprint( simplify(1+1/de)):
print(``):
print(`Proof: Consider  the integral`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+b*x)^(n+1),x=0..1)):
print(``):
print(`Let c be the constant  `):
print(log(1+b/a)/b):

print(``):

print(`It is readily seen that E(n) can be written as `):

print(``):

print(`E(n)=A(n)+B(n)*c `):
print(``):
print(`for some sequences A(n) and B(n) of RATIONAL NUMBERS `):

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

print(`It follows from the Almkvist-Zeilberger algorithm that`):
print(``):
print(`E(n), and hence A(n) and B(n), all satisfy the following second-order linear recurrence with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(``):
print(`but of course with different initial conditions.`):
print(``):



ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a,b)^n), E(n)=OMEGA( C2(a,b)^n ) `):
print(``):
print(`where C1(a,b) is the larger root and C2(a,b) the smaller root of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence`):
print(`that indicial polynomial, in N, is`):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
gu:=[solve(ope,N)]:

gu1:=evalf(subs({a=1,b=2},gu)):

if gu1[1]<gu1[2] then
 gu:=[gu[2],gu[1]]:
fi:

print(`Solving this quadratic, we get that the larger root, C1(a,b) is`):
print(``):
print(gu[1]):
print(``):
print(`and the smaller root, C2(a,b),  is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format,they are respectively`):
print(``):
lprint(gu[1],gu[2]):
print(``):

print(`We now need the following divisibilty lemma that is left as an exercise to the reader. `):
print(``):
print(`Lemma: A1(n)=lcm(1..n)*b^(2*n)*A(n) and B1(n)=lcm(1..n)*b^(2*n)*B(n) are always integers.`):
print(``):
print(`Let E1(n)=lcm(1..n)*b^(2*n)*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA((b^2*e*C1(a,b))^n ), E(n)=OMEGA( (b^2*e*C2(a,b))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log(b^2*e*C2(a,b))/log(b^2*e*C1(a,b)) `):
print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`It is readily seen that this is positive when a>(b-1/e)^2/4`):
print(``):
print(`and an irrationality measure can be taken to be 1+1/de, that establishes the theorem. QED`):
print(``):

end:


#GuessPK(P,x,N): Given a polynomial P in x of degree 1 or 2 and a positive integer N (at least 40) guesses the
#constant such that the coefficients of int(x^n*(1-x)^n/P^(n+1),x=0..1) when multiplied by lcm(1...n)*K^n integers. Try:
#GuessPK(1+x,x,40);
#GuessPK(1+x+x^2,x,40);
GuessPK:=proc(P,x,N) local gu,i,lu1,lu2:
gu:=ZugSeq(P,x,N)[2]:

lu1:=[seq(gu[i][1],i=1..N)]:
lu2:=[seq(gu[i][2],i=1..N)]:

[GuessK(lu1),GuessK(lu2)]:

end:




#FindGoodP(L,x): ouputs the list of [a+b*x+c*x^2,K,delta] such that the sequence 
#int( x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant
#int( 1/(a+b*x+c*x^2),x=0..1) for all 1<=a,b,c <= L. Try:
#FindGoodP(2);
FindGoodP:=proc(L,x) local lu,vu,gu,a,b,c,K,P,a1,b1,c1,K1,de:

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

gu:=[]:

for a1 from 1 to L do
 for b1 from 1 to L do
  for c1 from 1 to L do

  if igcd(a1,b1,c1)=1 and not type(sqrt(b1^2-4*a1*c1),integer) then

   P:=a1+b1*x+c1*x^2:
   vu:=GuessPK(P,x,40):

   if (vu[1]<>FAIL and vu[2]<>FAIL and    vu[1][2]<1000 and vu[2][2]<1000 and vu[1][1]=vu[2][1] ) then
     K1:=vu[1][1]:

de:=subs({a=a1,b=b1,c=c1,K=K1},lu):

    if evalf(de,20)>0 then
     gu:=[op(gu), [P,int(1/P,x=0..1),K1,de,evalf(1+1/de,20)]]:
    fi:
   fi:
 fi:
  od:
 od:
od:

gu:
end:


#FindGoodPAT(L,x): ouputs the list of [a+c*x^2,K,delta] such that the sequence 
#int( x^n*(1-x)^n/(a+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant
#int( 1/(a+c*x^2),x=0..1) for all 1<=a,c <= L. Try:
#FindGoodPAT(2);
FindGoodPAT:=proc(L) local lu,vu,gu,a,b,c,K,P,a1,b1,c1,K1,de:

lu:=subs(b=0,Delta2S(a,b,c,K)):

gu:=[]:

for a1 from 1 to L do
  for c1 from 1 to L do

  if igcd(a1,c1)=1 then

   P:=a1+c1*x^2:
   vu:=GuessPK(P,x,40):

   if (vu[1]<>FAIL and vu[2]<>FAIL and    vu[1][2]<1000 and vu[2][2]<1000 and vu[1][1]=vu[2][1] ) then
     K1:=vu[1][1]:

de:=subs({a=a1,c=c1,K=K1},lu):

    if evalf(de,20)>0 then
     gu:=[op(gu), [P,int(1/P,x=0..1),K1,de,evalf(1+1/de,20)]]:
    fi:
   fi:
 fi:
  od:
 od:

gu:
end:





#FindGoodPg(L): ouputs the list of [a+b*x+c*x^2,K,delta] such that the sequence 
#int( x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant
#int( 1/(a+b*x+c*x^2),x=0..1) for all 1<=a<=L, c<>0, and -L<=b,c <= L. Try:
#FindGoodPg(2);
FindGoodPg:=proc(L) local lu,vu,gu,a,b,c,K,P,a1,b1,c1,K1,de,x:

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

gu:=[]:

for a1 from 1 to L do
 for b1 from -L to L do
  for c1 from -L to L do


  if c1<>0 and  igcd(a1,b1,c1)=1  then

   P:=a1+b1*x+c1*x^2:
   if minimize(P,x=0..1)>0 then
   vu:=GuessPK(P,x,40):

   if (vu[1]<>FAIL and vu[2]<>FAIL and    vu[1][2]<1000 and vu[2][2]<1000 and vu[1][1]=vu[2][1] ) then
     K1:=vu[1][1]:

  de:=evalf(subs({a=a1,b=b1,c=c1,K=K1},lu),20):
    if de>0 then
     gu:=[op(gu), [P,int(1/P,x=0..1),K1,de]]:
    fi:
   fi:
 fi:
 fi:

  od:
 od:
od:

gu:
end:





#PaperQ(L): An article about the Diophantine approximations of Int(1/(a+b*x+c*x^2),x=0..1) for 1<=a,b,c<=L, and
#singling out those that yield irrationality proofs. Try:
#PaperQ(3);

PaperQ:=proc(L) local gu,E,n,C,ope,a,b,c,N,x,gu1,K,lu,de,i:
print(``):
print(`Irrationality Measures to some constants of the form int( 1/(a+b*x+c*x^2),x=0..1) with 1<=a,b,c <=`, L):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
gu:={}:

print(`Consider the constant`):
print(``):
print(C= Int(1/(a+b*x+c*x^2),x=0..1)):
print(``):
print(`with a,b,c, integers. We are interested in the diophantine approximations induced by the sequence`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1)):
print(``):
print(`that can be expressed as `):
print(``):
print(`E(n)=A(n)+B(n)*C`):
print(``):
print(`for sequences of rational numbers A(n) and B(n) `):
print(``):
ope:=AZd(x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x,n,N)[1]:
print(``):
print(`Using the Almkvist-Zeilberger algorithm it follows that E(n), and hence A(n) and B(n) satisfy`):
print(`the linear recurrence equation with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but of course with different initial conditions.`):
print(``):
ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a,b,c)^n), E(n)=OMEGA( C2(a,b,c)^n ) `):
print(``):
print(`where C1(a,b,c) is the larger root and C2(a,b,c) the smaller root of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence.`):
print(` That indicial polynomial, in N,  happens to be`):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):

gu:=[solve(ope,N)]:

gu:=[abs(gu[1]), abs(gu[2])]:

gu1:=evalf(subs({a=1,b=2,c=3},gu)):


if gu1[1]<gu1[2] then
 gu:=[gu[2],gu[1]]:
fi:

print(`Setting it equal to 0, and solving the quadratic equation, and taking the absolute value, we get that the large one,  C1(a,b,c), is`):
print(``):
print(gu[1]):
print(``):
print(`and the smaller one, C2(a,b,c) is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format, they are respectively, `):
print(``):
lprint(gu[1],gu[2]):
print(``):


print(`Suppose that we discover that there exists a constant K (that is either an integer or a square-root of an integer) such that `):
print(``):
print(` A1(n)=lcm(1..n)*K^n*A(n) and B1(n)=lcm(1..n)*K^n*B(n) are all integers.`):
print(``):
print(`(if K is a square-root of an integer only do it for even n .)`):
print(``):
print(`Let E1(n)=lcm(1..n)*K^n*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA(K^n*e*C1(a,b,c))^n ), E(n)=OMEGA( (K^n*e*C2(a,b,c))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log(K*e*C2(a,b,c))/log(K*e*C1(a,b,c)) `):
print(``):
de:=Delta2S(a,b,c,K):
print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`So it is good news whenever this is positive.`):
lu:=FindGoodP(L,x):

print(``):
print(`By searching for quadratic polynomials P with positive coefficients <=`, L, `we found the following list of `, nops(lu), `lucky cases. `):
print(``):
print(`For each of the lists below`):

print(`The first entry is the  polynomial in x, P(x), a certain irreducible quadratic.`):
print(``):
print(`the second entry is the constant whose irrationality we are claiming , namely the integral of 1/P(x) from x=0 to x=1.`):
print(``):
print(`The third entry is the magic constant K, mentioned above such that multiplying by K^n*lcm(1...n) makes everything integers.  `):
print(``):
print(`the fourth is the exact value of delta,`):
print(``):
print(` the fifth, and  last, entry is the implied irrationality measure, i.e. 1+1/delta,  in floating-point.`):
print(``):

for i from 1 to nops(lu) do
 print(``):
print(lu[i]):
 print(``):
od:

print(`and in Maple format `):

for i from 1 to nops(lu) do
 print(``):
lprint(lu[i]):
 print(``):
od:




end:




#PaperAT(L): An article about the Diophantine approximations of Int(1/(a+c*x^2),x=0..1) for 1<=a,c<=L that are
#successful
#PaperAT(L);

PaperAT:=proc(L) local gu,E,n,C,ope,a,b,c,N,x,gu1,K,lu,de,i:
print(``):
print(`Irrationality Measures to some constants of the form int( 1/(a+c*x^2),x=0..1) with 1<=a,c <=`, L):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
gu:={}:

print(`Consider the constant`):
print(``):
print(C= Int(1/(a+c*x^2),x=0..1)):
print(``):
print(`with a,c, integers. We are interested in the diophantine approximations induced by the sequence`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+c*x^2)^(n+1),x=0..1)):
print(``):
print(`that can be expressed as `):
print(``):
print(`E(n)=A(n)+B(n)*C`):
print(``):
print(`for sequences of rational numbers A(n) and B(n) `):
print(``):
ope:=AZd(x^n*(1-x)^n/(a+c*x^2)^(n+1),x,n,N)[1]:
print(``):
print(`Using the Almkvist-Zeilberger algorithm it follows that E(n), and hence A(n) and B(n) satisfy`):
print(`the linear recurrence equation with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but of course with different initial conditions.`):
print(``):
ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a,c)^n), E(n)=OMEGA( C2(a,c)^n ) `):
print(``):
print(`where C1(a,c) is the larger root and C2(a,c) the smaller root of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence.`):
print(` That indicial polynomial, in N,  happens to be`):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
gu:=[solve(ope,N)]:

gu1:=evalf(subs({a=1,c=3},gu)):

if gu1[1]<gu1[2] then
 gu:=[gu[2],gu[1]]:
fi:

print(`Setting it equal to 0, and solving the quadratic equation, we get that the larger root, C1(a,c) is`):
print(``):
print(gu[1]):
print(``):
print(`and the smaller root, C2(a,c) is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format, they are respectively, `):
print(``):
lprint(gu[1],gu[2]):
print(``):


print(`Suppose that we discover that there exists a constant K (that is either an integer or a square-root of an integer) such that `):
print(``):
print(` A1(n)=lcm(1..n)*K^n*A(n) and B1(n)=lcm(1..n)*K^n*B(n) are all integers.`):
print(``):
print(`(if K is a square-root of an integer only do it for even n .)`):
print(``):
print(`Let E1(n)=lcm(1..n)*K^n*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA(K^n*e*C1(a,c))^n ), E(n)=OMEGA( (K^n*e*C2(a,c))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log(K*e*C2(a,c))/log(K*e*C1(a,c)) `):
print(``):
de:=subs(b=0,Delta2S(a,b,c,K)):
print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`So it is good news whenever this is positive.`):
lu:=FindGoodPAT(L,x):

print(``):
print(`By searching for quadratic polynomials P with positive coefficients <=`, L, `we found the following list of `, nops(lu), `lucky cases. `):
print(``):
print(`For each of the lists below`):

print(`The first entry is the  polynomial in x, P(x), a certain irreducible quadratic.`):
print(``):
print(`the second entry is the constant whose irrationality we are claiming , namely the integral of 1/P(x) from x=0 to x=1.`):
print(``):
print(`The third entry is the magic constant K, mentioned above such that multiplying by K^n*lcm(1...n) makes everything integers.  `):
print(``):
print(`the fourth is the exact value of delta,`):
print(``):
print(` the fifth, and  last, entry is the implied irrationality measure, i.e. 1+1/delta,  in floating-point.`):
print(``):

for i from 1 to nops(lu) do
 print(``):
print(lu[i]):
 print(``):
od:

print(`and in Maple format `):

for i from 1 to nops(lu) do
 print(``):
lprint(lu[i]):
 print(``):
od:




end:



#Check37(L): checks that the conjecture that for all integers a that are 3,7 mod 8 up to L
#the divisibility constant of x^2+a is 2*sqrt(a). It outputs the exceptios. It the output
#is the empty set it is good news. Try:
#Check37(10);
Check37:=proc(K) local a,gu,x:

gu:={}:

for a from 3 to K do
 if member(a mod 8,{3,7}) then
   if GuessPK(x^2+a,x,60)[1][1]<>2*sqrt(a) then
     gu:=gu union {a}:
   fi:
 fi:
od:
gu:

end:



#Tovim(L): all the a from 1 to L such that
#the divisibility constant of x^2+a is 2*sqrt(a)

Tovim:=proc(L) local a,gu,x:

gu:={}:

for a from 3 to L do
   if GuessPK(x^2+a,x,60)[1][1]=2*sqrt(a) then
     gu:=gu union {a}:
   fi:
od:
gu:

end:

   










#PaperATS(): A computer-generated paper about the irrationality of 
#arctan(sqrt(a))/ sqrt(a), alias, int(1/(x^2+a),x=0..1) for ALL positive integers a that are 3 mod 8 or 7 mod 8.
#It gives a symbolic expression in a for the irrationality measure, complete with proof
#(modulo a divisibility lemma left to the reader). Try:
#PaperATS();
PaperATS:=proc() local E,n,a,b,c,de,ope,N,i,x,gu,gu1,lu:
print(`Irrationality Proofs and Measures for`, arctan(sqrt(a))/sqrt(a), ` for ALL positive integers a that are `):
print(`Congruent to 3  Modulo 4 `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

de:=subs({b=0,c=1},Delta2S(a,b,c,2*sqrt(a))):

print(`Theorem: Let a  be a positive integer such that a (mod 8)=3 or  a(mod 8)=7, then   arctan(sqrt(a))/sqrt(a)  is an irrational number`):
print(`and an irrationality measure is`, simplify(1+1/de) ):
print(``):
print(`In Maple format the irrationality measure is`):
print(``):
lprint( simplify(1+1/de)):
print(``):

print(`Proof: Consider  the integral`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+x^2)^(n+1),x=0..1)):
print(``):
print(`Let c be the constant  `):
print(arctan(sqrt(a))/sqrt(a)):

print(``):

print(`It is readily seen that E(n) can be written as `):

print(``):

print(`E(n)=A(n)+B(n)*c `):
print(``):
print(`for some sequences A(n) and B(n) of RATIONAL NUMBERS. `):


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


print(`It follows from the Almkvist-Zeilberger algorithm that`):
print(``):
print(`E(n), and hence A(n) and B(n), all satisfy the following second-order linear recurrence with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(``):
print(`but of course with different initial conditions.`):
print(``):



ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a)^n), E(n)=OMEGA( C2(a)^n ) `):
print(``):

print(`where C1(a) is the larger root and C2(a) the smaller root (in absolute value) of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence`):
print(`that indicial polynomial, in N, is`):
print(``):
print(ope):

print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
gu:=[solve(ope,N)]:

gu1:=evalf(subs({a=7},gu)):

if abs(gu1[1])<abs(gu1[2]) then
 gu:=[gu[2],gu[1]]:
fi:

print(`Solving this quadratic, we get that the absolute value of the root that has the larger absolute value, let's call is C1(a) is`):
print(``):
print(abs(gu[1])):
print(``):
print(`and the smaller root (in absolute value), let's call it, C2(a),  is`):
print(``):
print(abs(gu[2])):
print(``):
print(`and in Maple format,they are respectively`):
print(``):
lprint(-gu[1],gu[2]):
print(``):

print(`We now need the following divisibility lemma that is left as an exercise to the reader. `):
print(``):
print(`Lemma: A1(n)=lcm(1..n)*2^n*a^(trunc(n/2))*A(n) and  B1(n)=lcm(1..n)*2^n*a^(trunc(n/2))*B(n)   are always integers.`):
print(``):
print(`Let E1(n)=lcm(1..n)*2^n*a^(trunc(n/2))*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA((2*sqrt(a)*e*C1(a))^n ), E(n)=OMEGA( (2*sqrt(a)*e*C2(a))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log((2*sqrt(a)*e*C2(a))/log(b^2*e*C1(a,b)) `):

print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`It is readily seen that this is always positive`):
print(``):
lu:=simplify(1+1/de):
print(`and an irrationality measure can be taken to be 1+1/de that equals`):
print(``):
print(lu):
print(``):
print(`as claimed in the statement. `):
print(``):
print(`Finally, here is a plot ot this irrationality measure from a=3 to 1000.`):
print(``):
print(`Recall that it is ONLY valid for integers a that are 3(mod 8) or 7 (mod 8) `):
print(``):
print(plot(lu,a=3..1000)):
print(``):

print(``):

end:


#Bas(P,x): inputs a cubic polynomial in x with integer coefficients and outputs two constants {theta1,theta2} such that
#APx(F,x,n1) is always a linear combination of {1,theta1,theta2} with rational coefficients. Try:
#Bas(1+x+x^2+x^3,x);
Bas:=proc(P,x) local lu1,lu2,i:

if not (degree(P,x)=3 and {seq(type(coeff(P,x,i),integer),i=0..3)}={true})  then
  print(P, `must be a cubic polynomial in `, x, `with integer coefficients `):
 RETURN(FAIL):
fi:

lu1:=APx(P,x,1): lu1:=lu1-Frat(lu1):

lu2:=APx(P,x,2): lu2:=lu2-Frat(lu2):

{lu1,lu2}:


end:



 
#APxFS(P,x,n1,X): Like APxF(P,x,n1) but in terms of the initial conditions X[0], ..., X[ORD-1] where ORD is the order of the recurrence
#Try:
#APxFS(1+x+x^2+x^3,x,10,X);
APxFS:=proc(P,x,n1,X) local ope,N,ORD1,n,i1:
option remember:
ope:=AZd1(x^n*(1-x)^n/P^(n+1),x,n,N):
ORD1:=-ldegree(ope,N):

if n1<=ORD1-1 then
X[n1]:
else
add(subs(n=n1,coeff(ope,N,-i1))*APxFS(P,x,n1-i1,X),i1=1..ORD1):
fi:

end:


#ZugSeq(P,x,N): given a polynomial P in x of degree 1 or 2, output sAPx(P,x,0) (the irrational number we wish to approximate), let't call it c
#and the sequence of pairs [An,Bn] from n=1 to N such that APx(P,x,n)=An+Bn*c. Try:
#ZugSeq(1+x,x,30);
ZugSeq:=proc(P,x,N) local c,n1,lu,mu1,mu2,gu:

if not (degree(P,x)=1 or degree(P,x)=2) then
 print(P, `should be a polynomial in`, x, `of degree 1 or 2 `):
 RETURN(FAIL):
fi:
c:=APx(P,x,0):
gu:=[]:

for n1 from 1 to N do
 lu:=APxF(P,x,n1):
 mu1:=Frat(lu):
 mu2:=simplify((lu-mu1)/c):

 if not (type(mu2,fraction) or type(mu2,integer)) then
  RETURN(FAIL):
 fi:

 if simplify(lu-mu1-mu2*c)<>0 then
    RETURN(FAIL):
fi:

 gu:=[op(gu),[mu1,mu2]]:
od:

[c,gu]:

end:



#AppSeq(P,x,N): given a polynomial P in x of degree 1 or 2, output c= APx(P,x,0) (the irrational number we wish to approximate), 
#and the sequence of diophantine approximations inspired by APx(P,x,n); Try:
#AppSeq,30);
AppSeq:=proc(P,x,N) local gu,c, i:

gu:=ZugSeq(P,x,N):

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

c:=gu[1]:

gu:=gu[2]:

gu:=[seq(-gu[i][1]/gu[i][2],i=1..N)]:
[c,gu]:
end:


#The sequence of deltas for AppSeq from N-20 to N. Try:
#deltP(1+x,x,50);
deltP:=proc(P,x,N) local gu,c,i:


gu:=AppSeq(P,x,N):

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

if N<=20 then
 print(N, `must be at least 20 `):
fi:

c:=gu[1]:

gu:=gu[2]:

[seq(deltE(gu[i],c),i=N-20..N)]:

end:

#Hope(P,x,N): inputs a polynomial P of x with integer coefficients of degree 1 or 2 in x and a positive integer
#outputs the empircal delta for the N-th approximation. If it is positive there is hope that 
#we have an irratioality proof for the given constant. 
#The output is the constant and the empircal delta. Try:
#Hope(1+x,x,100);
Hope:=proc(P,x,N) local gu,c:

gu:=AppSeq(P,x,N):

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

c:=gu[1]:
gu:=gu[2]:

[c,deltE(gu[N],c)]:

end:

#GuessK1(Lis): Inputs a list of rational numbers, Lis, and outputs a positive integer K, and a small integer L such that
#Lis[n]*lcm(1...n)*K^n*L are all integers. Try
#GuessK1([seq(5^i/(dn(i)*7), i=1..30)]);
GuessK1:=proc(Lis) local N,n1,gu,K,i,mu,vu:

N:=nops(Lis):
if N<=21 then
 print(`The list has to be at least of length 21`):
fi:

gu:=ifactors(denom(Lis[N]*dn(N)))[2]:


mu:=1:

for i from 1 to nops(gu) do
 vu:=round(gu[i][2]/N):

  if abs(vu-gu[i][2]/N)<1/2 then
   mu:=mu*gu[i][1]^vu:
  fi:
od:
 

K:=mu:


gu:=lcm(seq(denom(Lis[n1]*dn(n1)*K^n1),n1=20..N)):


if {seq(denom(gu*Lis[n1]*dn(n1)*K^n1),n1=20..N)}<>{1} then
 RETURN(FAIL):
fi:

[K, gu]:


end:



#GuessK(Lis): Inputs a list of rational numbers, Lis, and outputs a positive integer K, and a small integer L such that
#Lis[n]*lcm(1...n)*K^n*L are all integers. Try
#GuessK([seq(5^i/(dn(i)*7), i=1..30)]);
GuessK:=proc(Lis) local N,K,gu,i,lu,ru:
K:=GuessK1(Lis):

if K[2]<=1000 then

if max(seq(denom(K[1]^i*Lis[i]*dn(i)*K[2]),i=1..nops(Lis)))>1000 then
  RETURN(FAIL):
 else
 RETURN(K):
fi:

fi:

N:=nops(Lis):

gu:=ifactors(denom(Lis[N]*dn(N)*K[1]^N))[2]:



lu:=1:

for i from 1 to nops(gu) do
 if abs(gu[i][2]/N-1/2)<1/5 then
  lu:=lu*gu[i][1]:
 fi:
od:


ru:=max(seq(denom(Lis[i]*dn(i)*K[1]^i*lu^(trunc(i/2))),i=1..nops(Lis))):



[K[1]*sqrt(lu),ru]:
end:








#deltT(P,x): inputs a polynomial P of degree 1 or 2, finds the theoretical delta for   the constant int(1/P,x=0..1) implied by
#the corresponding sequence int(x^n*(1-x)^n/P^(n+1),x=0..1); It returns, c, the conjectured K, and the delta, and its floating-point.
#So the output is
#[c,K,delta,evalf(delta), ), SmallestEmpricalDelta(from 180 to 200) ]: Try:
#deltT(1+x,x);
deltT:=proc(P,x) local gu,K1,K2,K,ope,c,N,i,n,lu:

gu:=ZugSeq(P,x,100):
c:=gu[1]:
gu:=gu[2]:



K1:=GuessK([seq(gu[i][1]*dn(i),i=1..nops(gu))]):

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

K2:=GuessK([seq(gu[i][2]*dn(i),i=1..nops(gu))]):

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




if K1[2]>1000 or K2[2]>1000 then
 RETURN(FAIL):
fi:

if K1[1]<>K2[1] then
 print(`The K-s are not the same`):
fi:

if type(K1[1],integer) and type(K2[1],integer) then
 K:=lcm(K1[1],K2[1]):
elif type(K1[1]^2,integer) and type(K2[1]^2,integer) then
 K:=sqrt(lcm(K1[1]^2,K2[1]^2)):
else
 RETURN(FAIL):
fi:




ope:=AZd(x^n*(1-x)^n/P^(n+1),x,n,N)[1]:




if degree(ope,N)<>2 then
 print(`The operator`, ope, `is not second-order `):
 RETURN(FAIL):
fi:

ope:=lcoeff(ope,n):



gu:=[solve(ope,N)]:

gu:=[abs(gu[1]), abs(gu[2])]:

if evalf(gu[1])>evalf(gu[2]) then
 gu:=[gu[2],gu[1]]:
fi:


lu:=-(log(gu[1])+log(K)+1)/(log(gu[2])+log(K)+1):

[c,K,lu,evalf(lu,20), min(op(deltP(P,x,100))) ]:

end:





#Delta1S(a,b):The explicit expression, in terms of a and b,  for the delta in the approximation of
#log((a+b)/a) implied by the integral int(x^n*(1-x)^n)/(a+b*x)^(n+1),x=0..1). Try:
#Delta1S(a,b);
Delta1S:=proc(a,b) local ope,gu,gu1,lu,x,n,N:

ope:=AZd((x^n*(1-x)^n)/(a+b*x)^(n+1),x,n,N)[1];

ope:=lcoeff(ope,n):

gu:=[solve(ope,N)]:


gu1:=subs({a=2,b=3},gu):

if evalf(abs(gu1[1]))>evalf(abs(gu1[2])) then
 gu:=[gu[2],gu[1]]:
fi:

lu:=-(log(gu[1])+2*log(b)+1)/(log(gu[2])+2*log(b)+1):

lu:
end:




#Delta2S(a,b,c,K):The explicit expression, in terms of a and b,  for the delta in the approximation of
# int(1/(a+b*x+c*x^2),x=0..1), implied by the integral int(x^n*(1-x)^n)/(a+b*x+c*x^2)^(n+1),x=0..1). 
#assuming that the divisibility constant is K. This is SYMBOLIC. To implement it you need to know K.
#Delta2S(a,b,c,K);
Delta2S:=proc(a,b,c,K) local ope,gu,gu1,lu,x,n,N:

ope:=AZd((x^n*(1-x)^n)/(a+b*x+c*x^2)^(n+1),x,n,N)[1];

ope:=lcoeff(ope,n):

gu:=[solve(ope,N)]:

gu:=[abs(gu[1]),abs(gu[2])]:

gu1:=subs({a=2,b=3,c=4},gu):

if evalf(gu1[1])>evalf(gu1[2]) then
 gu:=[gu[2],gu[1]]:
fi:

lu:=-(log(gu[1])+log(K)+1)/(log(gu[2])+log(K)+1):

lu:
end:


#CutO(a1): Inputs a and outputs the largest b such that Delta1S(a,b) is positive. Try:
#CutO(5);
CutO:=proc(a1) local lu,a,b,b1:
 lu:=Delta1S(a,b):

  for b1 from 1 while evalf(subs({a=a1,b=b1},lu))>0 do od:
b1-1:
end:

#CutOa(b): the smallest a such that Delta1S(a,b) is positive. Try:
#CutOa(b);
CutOa:=proc(b): ((b-exp(-1))/2)^2:end:

#CutOb(a): the largest b such that Delta1S(a,b) is positive. Try:
#CutOb(a);
CutOb:=proc(a) : 
(2*sqrt(a*exp(1))+1)/exp(1):
end:


#PaperA(): A computer-generated paper about the irrationality of log(1+b/a) for integer b< (2*sqrt(a*e)+1)/e
#and an explicit irratinality measure. It is a computer-generated version of Theorem 1 in  
#K. Alladi and M.L.Robinson, Legendre polynomials and irrationality, Crelle J.  318 (1980), 134-155.
#but is more stream-lined, direct and does NOT mention  Legedre polynomials. Type:
#PaperA();
PaperA:=proc() local E,n,a,b,de,ope,N,i,x,gu,gu1:
print(`Irrationality Proofs and Measures for log(1+b/a) for positive integers a and b such that a>(b-1/e)^2/4 `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

de:=Delta1S(a,b):

print(`Theorem: Let a and b be positive integers such that a>(b-1/e)^2/4, then log(1+b/a) is an irrational number`):
print(`and an irrationality measure is`, simplify(1+1/de) ):
print(``):
print(`In Maple format the irrationality measure is`):
print(``):
lprint( simplify(1+1/de)):
print(``):
print(`Proof: Consider  the integral`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+b*x)^(n+1),x=0..1)):
print(``):
print(`Let c be the constant  `):
print(log(1+b/a)/b):

print(``):

print(`It is readily seen that E(n) can be written as `):

print(``):

print(`E(n)=A(n)+B(n)*c `):
print(``):
print(`for some sequences A(n) and B(n) of RATIONAL NUMBERS `):

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

print(`It follows from the Almkvist-Zeilberger algorithm that`):
print(``):
print(`E(n), and hence A(n) and B(n), all satisfy the following second-order linear recurrence with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(``):
print(`but of course with different initial conditions.`):
print(``):



ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a,b)^n), E(n)=OMEGA( C2(a,b)^n ) `):
print(``):
print(`where C1(a,b) is the larger root and C2(a,b) the smaller root of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence`):
print(`that indicial polynomial, in N, is`):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
gu:=[solve(ope,N)]:

gu1:=evalf(subs({a=1,b=2},gu)):

if gu1[1]<gu1[2] then
 gu:=[gu[2],gu[1]]:
fi:

print(`Solving this quadratic, we get that the larger root, C1(a,b) is`):
print(``):
print(gu[1]):
print(``):
print(`and the smaller root, C2(a,b),  is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format,they are respectively`):
print(``):
lprint(gu[1],gu[2]):
print(``):

print(`We now need the following divisibilty lemma that is left as an exercise to the reader. `):
print(``):
print(`Lemma: A1(n)=lcm(1..n)*b^(2*n)*A(n) and B1(n)=lcm(1..n)*b^(2*n)*B(n) are always integers.`):
print(``):
print(`Let E1(n)=lcm(1..n)*b^(2*n)*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA((b^2*e*C1(a,b))^n ), E(n)=OMEGA( (b^2*e*C2(a,b))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log(b^2*e*C2(a,b))/log(b^2*e*C1(a,b)) `):
print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`It is readily seen that this is positive when a>(b-1/e)^2/4`):
print(``):
print(`and an irrationality measure can be taken to be 1+1/de, that establishes the theorem. QED`):
print(``):

end:


#GuessPK(P,x,N): Given a polynomial P in x of degree 1 or 2 and a positive integer N (at least 40) guesses the
#constant such that the coefficients of int(x^n*(1-x)^n/P^(n+1),x=0..1) when multiplied by lcm(1...n)*K^n integers. Try:
#GuessPK(1+x,x,40);
#GuessPK(1+x+x^2,x,40);
GuessPK:=proc(P,x,N) local gu,i,lu1,lu2:
gu:=ZugSeq(P,x,N)[2]:

lu1:=[seq(gu[i][1],i=1..N)]:
lu2:=[seq(gu[i][2],i=1..N)]:

[GuessK(lu1),GuessK(lu2)]:

end:




#FindGoodP(L,x): ouputs the list of [a+b*x+c*x^2,K,delta] such that the sequence 
#int( x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant
#int( 1/(a+b*x+c*x^2),x=0..1) for all 1<=a,b,c <= L. Try:
#FindGoodP(2);
FindGoodP:=proc(L,x) local lu,vu,gu,a,b,c,K,P,a1,b1,c1,K1,de:

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

gu:=[]:

for a1 from 1 to L do
 for b1 from 1 to L do
  for c1 from 1 to L do

  if igcd(a1,b1,c1)=1 and not type(sqrt(b1^2-4*a1*c1),integer) then

   P:=a1+b1*x+c1*x^2:
   vu:=GuessPK(P,x,40):

   if (vu[1]<>FAIL and vu[2]<>FAIL and    vu[1][2]<1000 and vu[2][2]<1000 and vu[1][1]=vu[2][1] ) then
     K1:=vu[1][1]:

de:=subs({a=a1,b=b1,c=c1,K=K1},lu):

    if evalf(de,20)>0 then
     gu:=[op(gu), [P,int(1/P,x=0..1),K1,de,evalf(1+1/de,20)]]:
    fi:
   fi:
 fi:
  od:
 od:
od:

gu:
end:


#FindGoodPAT(L,x): ouputs the list of [a+c*x^2,K,delta] such that the sequence 
#int( x^n*(1-x)^n/(a+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant
#int( 1/(a+c*x^2),x=0..1) for all 1<=a,c <= L. Try:
#FindGoodPAT(2);
FindGoodPAT:=proc(L) local lu,vu,gu,a,b,c,K,P,a1,b1,c1,K1,de:

lu:=subs(b=0,Delta2S(a,b,c,K)):

gu:=[]:

for a1 from 1 to L do
  for c1 from 1 to L do

  if igcd(a1,c1)=1 then

   P:=a1+c1*x^2:
   vu:=GuessPK(P,x,40):

   if (vu[1]<>FAIL and vu[2]<>FAIL and    vu[1][2]<1000 and vu[2][2]<1000 and vu[1][1]=vu[2][1] ) then
     K1:=vu[1][1]:

de:=subs({a=a1,c=c1,K=K1},lu):

    if evalf(de,20)>0 then
     gu:=[op(gu), [P,int(1/P,x=0..1),K1,de,evalf(1+1/de,20)]]:
    fi:
   fi:
 fi:
  od:
 od:

gu:
end:





#FindGoodPg(L): ouputs the list of [a+b*x+c*x^2,K,delta] such that the sequence 
#int( x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1) yields an irrationality measure for the constant
#int( 1/(a+b*x+c*x^2),x=0..1) for all 1<=a<=L, c<>0, and -L<=b,c <= L. Try:
#FindGoodPg(2);
FindGoodPg:=proc(L) local lu,vu,gu,a,b,c,K,P,a1,b1,c1,K1,de,x:

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

gu:=[]:

for a1 from 1 to L do
 for b1 from -L to L do
  for c1 from -L to L do


  if c1<>0 and  igcd(a1,b1,c1)=1  then

   P:=a1+b1*x+c1*x^2:
   if minimize(P,x=0..1)>0 then
   vu:=GuessPK(P,x,40):

   if (vu[1]<>FAIL and vu[2]<>FAIL and    vu[1][2]<1000 and vu[2][2]<1000 and vu[1][1]=vu[2][1] ) then
     K1:=vu[1][1]:

  de:=evalf(subs({a=a1,b=b1,c=c1,K=K1},lu),20):
    if de>0 then
     gu:=[op(gu), [P,int(1/P,x=0..1),K1,de]]:
    fi:
   fi:
 fi:
 fi:

  od:
 od:
od:

gu:
end:





#PaperQ(L): An article about the Diophantine approximations of Int(1/(a+b*x+c*x^2),x=0..1) for 1<=a,b,c<=L, and
#singling out those that yield irrationality proofs. Try:
#PaperQ(3);

PaperQ:=proc(L) local gu,E,n,C,ope,a,b,c,N,x,gu1,K,lu,de,i:
print(``):
print(`Irrationality Measures to some constants of the form int( 1/(a+b*x+c*x^2),x=0..1) with 1<=a,b,c <=`, L):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
gu:={}:

print(`Consider the constant`):
print(``):
print(C= Int(1/(a+b*x+c*x^2),x=0..1)):
print(``):
print(`with a,b,c, integers. We are interested in the diophantine approximations induced by the sequence`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x=0..1)):
print(``):
print(`that can be expressed as `):
print(``):
print(`E(n)=A(n)+B(n)*C`):
print(``):
print(`for sequences of rational numbers A(n) and B(n) `):
print(``):
ope:=AZd(x^n*(1-x)^n/(a+b*x+c*x^2)^(n+1),x,n,N)[1]:
print(``):
print(`Using the Almkvist-Zeilberger algorithm it follows that E(n), and hence A(n) and B(n) satisfy`):
print(`the linear recurrence equation with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but of course with different initial conditions.`):
print(``):
ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a,b,c)^n), E(n)=OMEGA( C2(a,b,c)^n ) `):
print(``):
print(`where C1(a,b,c) is the larger root and C2(a,b,c) the smaller root of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence.`):
print(` That indicial polynomial, in N,  happens to be`):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):

gu:=[solve(ope,N)]:

gu:=[abs(gu[1]), abs(gu[2])]:

gu1:=evalf(subs({a=1,b=2,c=3},gu)):


if gu1[1]<gu1[2] then
 gu:=[gu[2],gu[1]]:
fi:

print(`Setting it equal to 0, and solving the quadratic equation, and taking the absolute value, we get that the large one,  C1(a,b,c), is`):
print(``):
print(gu[1]):
print(``):
print(`and the smaller one, C2(a,b,c) is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format, they are respectively, `):
print(``):
lprint(gu[1],gu[2]):
print(``):


print(`Suppose that we discover that there exists a constant K (that is either an integer or a square-root of an integer) such that `):
print(``):
print(` A1(n)=lcm(1..n)*K^n*A(n) and B1(n)=lcm(1..n)*K^n*B(n) are all integers.`):
print(``):
print(`(if K is a square-root of an integer only do it for even n .)`):
print(``):
print(`Let E1(n)=lcm(1..n)*K^n*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA(K^n*e*C1(a,b,c))^n ), E(n)=OMEGA( (K^n*e*C2(a,b,c))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log(K*e*C2(a,b,c))/log(K*e*C1(a,b,c)) `):
print(``):
de:=Delta2S(a,b,c,K):
print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`So it is good news whenever this is positive.`):
lu:=FindGoodP(L,x):

print(``):
print(`By searching for quadratic polynomials P with positive coefficients <=`, L, `we found the following list of `, nops(lu), `lucky cases. `):
print(``):
print(`For each of the lists below`):

print(`The first entry is the  polynomial in x, P(x), a certain irreducible quadratic.`):
print(``):
print(`the second entry is the constant whose irrationality we are claiming , namely the integral of 1/P(x) from x=0 to x=1.`):
print(``):
print(`The third entry is the magic constant K, mentioned above such that multiplying by K^n*lcm(1...n) makes everything integers.  `):
print(``):
print(`the fourth is the exact value of delta,`):
print(``):
print(` the fifth, and  last, entry is the implied irrationality measure, i.e. 1+1/delta,  in floating-point.`):
print(``):

for i from 1 to nops(lu) do
 print(``):
print(lu[i]):
 print(``):
od:

print(`and in Maple format `):

for i from 1 to nops(lu) do
 print(``):
lprint(lu[i]):
 print(``):
od:




end:




#PaperAT(L): An article about the Diophantine approximations of Int(1/(a+c*x^2),x=0..1) for 1<=a,c<=L that are
#successful
#PaperAT(L);

PaperAT:=proc(L) local gu,E,n,C,ope,a,b,c,N,x,gu1,K,lu,de,i:
print(``):
print(`Irrationality Measures to some constants of the form int( 1/(a+c*x^2),x=0..1) with 1<=a,c <=`, L):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
gu:={}:

print(`Consider the constant`):
print(``):
print(C= Int(1/(a+c*x^2),x=0..1)):
print(``):
print(`with a,c, integers. We are interested in the diophantine approximations induced by the sequence`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+c*x^2)^(n+1),x=0..1)):
print(``):
print(`that can be expressed as `):
print(``):
print(`E(n)=A(n)+B(n)*C`):
print(``):
print(`for sequences of rational numbers A(n) and B(n) `):
print(``):
ope:=AZd(x^n*(1-x)^n/(a+c*x^2)^(n+1),x,n,N)[1]:
print(``):
print(`Using the Almkvist-Zeilberger algorithm it follows that E(n), and hence A(n) and B(n) satisfy`):
print(`the linear recurrence equation with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(`but of course with different initial conditions.`):
print(``):
ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a,c)^n), E(n)=OMEGA( C2(a,c)^n ) `):
print(``):
print(`where C1(a,c) is the larger root and C2(a,c) the smaller root of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence.`):
print(` That indicial polynomial, in N,  happens to be`):
print(``):
print(ope):
print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
gu:=[solve(ope,N)]:

gu1:=evalf(subs({a=1,c=3},gu)):

if gu1[1]<gu1[2] then
 gu:=[gu[2],gu[1]]:
fi:

print(`Setting it equal to 0, and solving the quadratic equation, we get that the larger root, C1(a,c) is`):
print(``):
print(gu[1]):
print(``):
print(`and the smaller root, C2(a,c) is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format, they are respectively, `):
print(``):
lprint(gu[1],gu[2]):
print(``):


print(`Suppose that we discover that there exists a constant K (that is either an integer or a square-root of an integer) such that `):
print(``):
print(` A1(n)=lcm(1..n)*K^n*A(n) and B1(n)=lcm(1..n)*K^n*B(n) are all integers.`):
print(``):
print(`(if K is a square-root of an integer only do it for even n .)`):
print(``):
print(`Let E1(n)=lcm(1..n)*K^n*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA(K^n*e*C1(a,c))^n ), E(n)=OMEGA( (K^n*e*C2(a,c))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log(K*e*C2(a,c))/log(K*e*C1(a,c)) `):
print(``):
de:=subs(b=0,Delta2S(a,b,c,K)):
print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`So it is good news whenever this is positive.`):
lu:=FindGoodPAT(L,x):

print(``):
print(`By searching for quadratic polynomials P with positive coefficients <=`, L, `we found the following list of `, nops(lu), `lucky cases. `):
print(``):
print(`For each of the lists below`):

print(`The first entry is the  polynomial in x, P(x), a certain irreducible quadratic.`):
print(``):
print(`the second entry is the constant whose irrationality we are claiming , namely the integral of 1/P(x) from x=0 to x=1.`):
print(``):
print(`The third entry is the magic constant K, mentioned above such that multiplying by K^n*lcm(1...n) makes everything integers.  `):
print(``):
print(`the fourth is the exact value of delta,`):
print(``):
print(` the fifth, and  last, entry is the implied irrationality measure, i.e. 1+1/delta,  in floating-point.`):
print(``):

for i from 1 to nops(lu) do
 print(``):
print(lu[i]):
 print(``):
od:

print(`and in Maple format `):

for i from 1 to nops(lu) do
 print(``):
lprint(lu[i]):
 print(``):
od:




end:



#Check37(L): checks that the conjecture that for all integers a that are 3,7 mod 8 up to L
#the divisibility constant of x^2+a is 2*sqrt(a). It outputs the exceptios. It the output
#is the empty set it is good news. Try:
#Check37(10);
Check37:=proc(K) local a,gu,x:

gu:={}:

for a from 3 to K do
 if member(a mod 8,{3,7}) then
   if GuessPK(x^2+a,x,60)[1][1]<>2*sqrt(a) then
     gu:=gu union {a}:
   fi:
 fi:
od:
gu:

end:



#Tovim(L): all the a from 1 to L such that
#the divisibility constant of x^2+a is 2*sqrt(a)

Tovim:=proc(L) local a,gu,x:

gu:={}:

for a from 3 to L do
   if GuessPK(x^2+a,x,60)[1][1]=2*sqrt(a) then
     gu:=gu union {a}:
   fi:
od:
gu:

end:

   










#PaperATS(): A computer-generated paper about the irrationality of 
#arctan(sqrt(a))/ sqrt(a), alias, int(1/(x^2+a),x=0..1) for ALL positive integers a that are 3 mod 8 or 7 mod 8.
#It gives a symbolic expression in a for the irrationality measure, complete with proof
#(modulo a divisibility lemma left to the reader). Try:
#PaperATS();
PaperATS:=proc() local E,n,a,b,c,de,ope,N,i,x,gu,gu1,lu:
print(`Irrationality Proofs and Measures for`, arctan(sqrt(a))/sqrt(a), ` for ALL positive integers a that are `):
print(`Congruent to 3 Modulo 4 `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

de:=subs({b=0,c=1},Delta2S(a,b,c,2*sqrt(a))):

print(`Theorem: Let a  be a positive integer such that a (mod 4)=3,  then   arctan(sqrt(a))/sqrt(a)  is an irrational number`):
print(`and an irrationality measure is`, simplify(1+1/de) ):
print(``):
print(`In Maple format the irrationality measure is`):
print(``):
lprint( simplify(1+1/de)):
print(``):

print(`Proof: Consider  the integral`):
print(``):
print(E(n)=Int(x^n*(1-x)^n/(a+x^2)^(n+1),x=0..1)):
print(``):
print(`Let c be the constant  `):
print(arctan(sqrt(a))/sqrt(a)):

print(``):

print(`It is readily seen that E(n) can be written as `):

print(``):

print(`E(n)=A(n)+B(n)*c `):
print(``):
print(`for some sequences A(n) and B(n) of RATIONAL NUMBERS. `):


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


print(`It follows from the Almkvist-Zeilberger algorithm that`):
print(``):
print(`E(n), and hence A(n) and B(n), all satisfy the following second-order linear recurrence with polynomial coefficients`):
print(``):
print(add(coeff(ope,N,i)*F(n+i),i=0..degree(ope,N))=0):
print(``):
print(``):
print(`but of course with different initial conditions.`):
print(``):



ope:=lcoeff(ope,n):
ope:=normal(ope/coeff(ope,N,0)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
print(``): 
print(`It follows from the Poincarre lemma that `):
print(``):
print(`A(n),B(n)=OMEGA(C1(a)^n), E(n)=OMEGA( C2(a)^n ) `):
print(``):

print(`where C1(a) is the larger root and C2(a) the smaller root (in absolute value) of the indicial equation of the constant-coefficient`):
print(`recurrence approximating the above recurrence`):
print(`that indicial polynomial, in N, is`):
print(``):
print(ope):

print(``):
print(`and in Maple format`):
print(``):
lprint(ope):
print(``):
gu:=[solve(ope,N)]:

gu1:=evalf(subs({a=7},gu)):

if abs(gu1[1])<abs(gu1[2]) then
 gu:=[gu[2],gu[1]]:
fi:

print(`Solving this quadratic, we get that the absolute value of the root that has the larger absolute value, let's call is C1(a) is`):
print(``):
print(abs(gu[1])):
print(``):
print(`and the smaller root (in absolute value), let's call it, C2(a),  is`):
print(``):
print(abs(gu[2])):
print(``):
print(`and in Maple format,they are respectively`):
print(``):
lprint(-gu[1],gu[2]):
print(``):

print(`We now need the following divisibility lemma that is left as an exercise to the reader. `):
print(``):
print(`Lemma: A1(n)=lcm(1..n)*2^n*a^(trunc(n/2))*A(n) and  B1(n)=lcm(1..n)*2^n*a^(trunc(n/2))*B(n)   are always integers.`):
print(``):
print(`Let E1(n)=lcm(1..n)*2^n*a^(trunc(n/2))*E(n) , then `):
print(``):
print(`E1(n)=A1(n)+B1(n)*c `):
print(``):
print(`where now A1(n) and B1(n) are INTEGERS.`):
print(``):
print(`Since lcm(1..n)=OMEGA(e^n), we have `):
print(``):
print(`A1(n),B1(n)=OMEGA((2*sqrt(a)*e*C1(a))^n ), E(n)=OMEGA( (2*sqrt(a)*e*C2(a))^n ) `):
print(``):
print(`Hence `):
print(`E(n)=OMEGA(1/max(A(n),B(n))^delta ), where delta is `):
print(``):
print(`-log((2*sqrt(a)*e*C2(a))/log(b^2*e*C1(a,b)) `):

print(``):
print(`that equals `, de):
print(``):
print(`and in Maple format it is`):
print(``):
lprint(de):
print(``):
print(`It is readily seen that this is always positive`):
print(``):
lu:=simplify(1+1/de):
print(`and an irrationality measure can be taken to be 1+1/de that equals`):
print(``):
print(lu):
print(``):
print(`as claimed in the statement. `):
print(``):
print(`Finally, here is a plot ot this irrationality measure from a=3 to 1000.`):
print(``):
print(`Recall that it is ONLY valid for integers a that are 3(mod 8) or 7 (mod 8) `):
print(``):
print(plot(lu,a=3..1000)):
print(``):

print(``):

end:


#Bas(P,x): inputs a cubic polynomial in x with integer coefficients and outputs two constants {theta1,theta2} such that
#APx(F,x,n1) is always a linear combination of {1,theta1,theta2} with rational coefficients. Try:
#Bas(1+x+x^2+x^3,x);
Bas:=proc(P,x) local lu1,lu2,i:

if not (degree(P,x)=3 and {seq(type(coeff(P,x,i),integer),i=0..3)}={true})  then
  print(P, `must be a cubic polynomial in `, x, `with integer coefficients `):
 RETURN(FAIL):
fi:

lu1:=APx(P,x,1): lu1:=lu1-Frat(lu1):

lu2:=APx(P,x,2): lu2:=lu2-Frat(lu2):

{lu1,lu2}:


end:



 
#APxFS(P,x,n1,X): Like APxF(P,x,n1) but in terms of the initial conditions X[0], ..., X[ORD-1] where ORD is the order of the recurrence
#Try:
#APxFS(1+x+x^2+x^3,x,10,X);
APxFS:=proc(P,x,n1,X) local ope,N,ORD1,n,i1:
option remember:
ope:=AZd1(x^n*(1-x)^n/P^(n+1),x,n,N):
ORD1:=-ldegree(ope,N):

if n1<=ORD1-1 then
X[n1]:
else
add(subs(n=n1,coeff(ope,N,-i1))*APxFS(P,x,n1-i1,X),i1=1..ORD1):
fi:

end:



#TupSeq(P,x,N): given a polynomial P in x  output the sequence of coefficients of X[0], .., X[Order-1] in
#APxFS(P,x,n) from n=1 to N. Try:
#TupSeq(1+x+x^2+x^3,x,30,X);
TupSeq:=proc(P,x,N) local X,n,S,n1,i1,ope,ORD1:
ope:=AZd1(x^n*(1-x)^n/P^(n+1),x,n,S):
ORD1:=-ldegree(ope,S):

[seq([seq(coeff(APxFS(P,x,n1,X),X[i1],1),n1=1..N)],i1=0..ORD1-1)]:

end:

#TripSeqSp(N): the sequence of coefficients of 1, log(2), and Pi for APx(1+x+x^2+x^3,x,n) for n from 1 to N. Try:
#TripSeqSp(30);
TripSeqSp:=proc(N) local gu,x,n1:
gu:=[seq(APxF(1+x+x^2+x^3,x,n1),n1=1..N)]:

[
[seq(coeff(coeff(gu[n1],log(2),0),Pi,0),n1=1..N)],
[seq(coeff(gu[n1],log(2),1),n1=1..N)],
[seq(coeff(gu[n1],Pi,1),n1=1..N)]
]:

end:


#DeltSp(): The exponent such that APxF(1+x+x^2+x^3,x) gives a simultaneous approximation to 1, log(2), Pi. Try:
#DeltSp();
DeltSp:= proc() local gu,n1,ope,ope1,x,n,N:
gu:=TripSeqSp(40):

print(
{seq(denom(gu[1][n1]*dn(n1)*2^n1),n1=1..40)},
{seq(denom(gu[2][n1]*2^n1),n1=1..40)},
{seq(denom(gu[3][n1]*4^n1),n1=1..40)}):


ope:=AZd(x^n*(1-x)^n/(1+x+x^2+x^3)^(n+1),x,n,N)[1];

ope1:=lcoeff(ope,n):


ope1:=subs(N=sqrt(N),ope1):

gu:=[solve(ope1,N)]:

gu:=[sqrt(abs(gu[1])),sqrt(abs(gu[2]))]:

-(  log(gu[1])+log(4)+1)/(log(gu[2])+log(4)+1):


end: