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

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



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

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

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

 print(`For a list of the procedures taken from the Heimonen, Matala-Aho, Vaannanen paper type ezra2();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 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: Book, Paper, Prop1 `):
 print(``):

else
ezra(args):
fi:

end:

ezra2:=proc()

if args=NULL then
 print(` The  procedures from the Heimonen et. al. paper are:  IrrM, kstar, lamh, muh `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Akna, Aka1, c1, c2, ChA, CheckDH, ConjK, ConjK1,  COV,  Deco1, dn, EvalCF, GrabHP, GuessDiv, ReRa, Rnus, RP `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: AknaF,  AZd1, Bkna, deltE, deltP,  Info, Seqs , Tikva`):
 print(` `):


elif nops([args])=1 and op(1,[args])=Aka1 then
print(`Aka1(k,a,X): Expressing Akna(k,1,a) in terms of X=Akna(k,0,a); Try:`):
print(`Aka1(3,4,X);`):


elif nops([args])=1 and op(1,[args])=Akna then
print(`Akna(k,n,a): int(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x=0..1); Done DIRECTLY`):
print(`Try: `):
print(`Akna(3,4,4);`):

elif nops([args])=1 and op(1,[args])=AknaF then
print(`AknaF(k,n1,a,X): Fast version of Akna(k,n1,a) using the Almkvist-Zeilberger, where X denotes Akna(k,0,a), i.e.`):
print(`int(1/(1+x^k/a),x=0..1);`):
print(`AknaF(k,n1,a,X): Fast version of Akna(k,n1,a) using the Almkvist-Zeilberger algorithm. Try:`):
print(`AknaF(3,10,3,X);`):

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])=Bkna then
print(`Bkna(k,n,a): The pair of rational numbers [c(n),d(n)] such that Akna(k,n,a)=c(n)+d(n)*Akna(k,0,a) of FAIL.`):
print(`It also returns the ratio that approximates Akna(k,0,a), the error and the empirical delta`):
print(`Try:`):
print(`Bkna(3,3,1);`):


elif nops([args])=1 and op(1,[args])=Book then
print(`Book(k, Maxa, G,N): Inputs a positive integer k, a positive integer Maxa, and positive integers G and N (parameters for guessing)`):
print(`outputs a book about irrationaility, and irrationality measures of int(1/(1+x^k/a),x=0..1) for a from 1 to Maxa`):
print(`Try:`):
print(`Book(1,20,200,200);`):

elif nops([args])=1 and op(1,[args])=c1 then
print(`c1(a): The constant c1(a) such that |Akna(k,n,a)|<=C/c1(a)^n. Try:`):
print(`c1(a); c1(5);`):

elif nops([args])=1 and op(1,[args])=c2 then
print(`c2(k,a): The constant c2(k,a) such that |Bkna(k,n,a)[1]|<=C*c2(a,k)^n and |Bkna(k,n,a)[2]|<=C*c2(a,k)^n`):
print(`a is a pos. integer and k is a pos. integer. Try: `):
print(`c2(3,4);`):

elif nops([args])=1 and op(1,[args])=ChA then
print(`ChA(k,a,N0): Checks procedure AknaF(k,n,a,X) for n from 0 to N0. Try:`):
print(`ChA(2,2,20);`):

elif nops([args])=1 and op(1,[args])=CheckDH then
print(`CheckDH(k,a,N): Checks the divisibility hypothesis for Seqs(k,a,N). Try:`):
print(`CheckDH(7,11,100);`):




elif nops([args])=1 and op(1,[args])=ConjK then
print(`ConjK(k,a,N,G): Given a positive integer k and a rational number a, and a fairly large positive integer N,`):
print(`outputs the rational number K with numerators and denominators <=G `):
print(`such that the the sequences of numerators and denominators in the diophantine approximations to int(1/(1+x^k.a),x=0..1) `):
print(`are such that if you multiply the n-th entry by lcm(1..k*n) and divide by K^n you always get an integer, up to the N-th term.`):
print(`Try: `):
print(`ConjK(2,11,40,30);`):

elif nops([args])=1 and op(1,[args])=ConjK1 then
print(`ConjK1(L,G): inputs a sequence of rational numbers, L, and a limit G, outputs a rational number with numerator`):
print(`and denominator <=G such that L[i]*K^i are all integers. Try:`):
print(`ConjK1([seq((2/3)^i,i=1..20)],20);`):


elif nops([args])=1 and op(1,[args])=COV then
print(`COV(L,K1,K2): inputs a list of rational numbers L, a positive integer K1, and a rational number K2, outputs the sequence L[n]*dn(K1*n)*K2^n1`):
print(`Try:`):
print(`COV(Seqs(3,8,20)[1],3,1/8);`):

elif nops([args])=1 and op(1,[args])=Deco11 then
print(`Deco11(a,n,r): Given a rational number a, find the smallest K1 such that the denominator of a*dn(n*K1) has <=r factors`):
print(`Try:`):
print(`Deco11(1/dn(30)/5^30,10,3);`):

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 EMPIRICAL delta such that |a-c|=1/denom(a)^(1+delta). `):
print(`For example, try deltaE(22/7,evalf(Pi)):`):

elif nops([args])=1 and op(1,[args])=deltP then
print(`deltP(k,a,K1,K2): let alpha:=int(1/(1+x^k/a),x=0..1). Write`):
print(`int(x^n*(1-x)^n/(1+(x/a)^k)^(n+1),x=0..1)=A(n)-B(n)*alpha`):
print(`where A(n) and B(n) are specific rational numbers.`):
print(`Assume that you proved that A'(n):=A(n)*lcm(1..K1*n)/K2^n,  B'(n):=B(n)*lcm(1..K1*n)/K2^n, `):
print(`are integers, it outputs the number delta such that`):
print(`|alpha-A'(n)/B'(n)|<=CONSTANT/B'(n)^(1+delta). If delta is >0, then this proves the irrationality of`):
print(`alpha, and the irrationality measure is 1+1/delta. Try:`):
print(`deltP(3,2,2,3);`):

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

elif nops([args])=1 and op(1,[args])=EvalCF then
print(`EvalCF(L): evaluates the continued fraction. Try:`):
print(`EvalCF([5,4,3]);`):

elif nops([args])=1 and op(1,[args])=GrabHP then
print(`GrabHP(n): inputs an integer n and outputs [m,r] such that n/m^r is an integer and r is as large as possible. Try:`):
print(`GrabHP(10^20*101);`):

elif nops([args])=1 and op(1,[args])=GuessDiv then
print(`GuessDiv(L): Given a sequence of integers L, guesses an integer m, and a rational number such that`):
print(`L[n]/m^(trunc(n*alpha)) are integers. Try:`):
print(`GuessDiv([seq(5^(3*i)*13,i=1..10)]);`):


elif nops([args])=1 and op(1,[args])=Info then
print(`Info(k,a,G,N): Inputs a pos. integer k, a rational number a, a fairly large (say 100) pos. integer G, and a fairly large`):
print(`pos. integer N (say 200), outputs a list [constant, K,deltP,evalf(deltP), deltE] , where constant is the exact value of`):
print(`constant of interest, int(1/(1+x^k/a),x=0..1),`):
print(`K is the rational number such that `):
print(`if A(n) B(n) are the sequences of rational number such that A(n)/B(n) is a diophantine approximation to`):
print(`int(1/(1+x^k/a),x=0..1), we have conjecturally that both A(n)*K^n, and B(n)*K^n are always integers (proved, so far for n<=N` ):
print(`where the numerator and denominator of K are <=G. deltP is the implied rigorous (modulo the divisibility lemma ` ):
print(`and deltE is the smallest empirical delta from N/2 to N. Try:`):
print(`Info(2,11,100,200);`):

elif nops([args])=1 and op(1,[args])=IrrM then
print(`IrrM(k,a): The irrationality measure according to HMV with beta=1 . Try: `):
print(`IrrM(2,11);`):

elif nops([args])=1 and op(1,[args])=kstar then
print(`kstar(k,a): inputs a pos. integer k and a rational mumber a= -s/r, outputs the denomiator of (s-r)/k `):
print(`and (s-r)/(k*prod(p, p|k). Try:`):
print(`kstar(2,-1/7);`):

elif nops([args])=1 and op(1,[args])=lamh then
print(`lamh(h): h/phi(h)*sum(1/i,i=1..h, gcd(i,h)=1). Try:`):
print(`lamh(20);`):

elif nops([args])=1 and op(1,[args])=muh then
print(`muh(h): the product of p^(1/(p-1)) over all prime divisors of h. Try`):
print(`muh(11);`):

elif nops([args])=1 and op(1,[args])=Paper then
print(`Paper(k,a,G,N): A paper about the irrationality of int(1/(1+x^k/a),x=0..1); G and N are parameters for estimating and guessing. Try:`):
print(`Paper(2,11,100,200):`):


elif nops([args])=1 and op(1,[args])=Prop1 then
print(`Prop1(k,a,G,N,co): A proposition numbered co, about the irrationality of int(1/(1+x^k/a),x=0..1); G and N are parameters for estimating and guessing. Try:`):
print(`Prop1(2,11,100,200,10):`):

elif nops([args])=1 and op(1,[args])=ReRa then
print(`ReRa(a,K): inputs a real number approximates it with a rational number. K is a big number indicationg where`):
print(`to truncate the continued fraction. It also returns the error.`):
print(`It returns FAIL if it is fails.Try: `):
print(`ReRa(101/89+10^(-100),10^10);`):

elif nops([args])=1 and op(1,[args])=Rnus then
print(`Rnus(N): The first N positive rational numbers, Try:`):
print(`Rnus(10);`):


elif nops([args])=1 and op(1,[args])=RP then
print(`RP(a): given a real number that has one part that is a fraction, call it f, reutrns the pair [f,a-f]. Try: `):
print(`RP(5/3+sqrt(2));`):

elif nops([args])=1 and op(1,[args])=Seqs then
print(`Seqs(k,a,N): The two sequence B(k,n,a)[1] and B(k,n,a)[2] for n from 1 to N. Try:`):
print(`Seqs(3,2,20);`):

elif nops([args])=1 and op(1,[args])=Tikva then
print(`Tikva(k,a,N1,N2): Given k and a and a  positive integers 20<=N1<N2  outputs the smallest`):
print(`empirical delta for N1 to N2 for the sequence of rational numbers approximating int(1/(1+x^k/a),x=0..1);`):
print(`Tikva(1,1,30,50);`):

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

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

 


##############start generalized arctan(1/a)


#Akna(k,n,a): int(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x=0..1);
#Try:
#Akna(k,n,a);
Akna:=proc(k,n,a) local x:
option remember:
int(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x=0..1);
end:

#RP(a): given a real number that has one part that is a fraction, call it f, reutrns the pair [f,a-f]. Try
#RP(5/3+sqrt(2));
RP:=proc(a) local i,f:

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

for i from 1 to nops(a) do
 if type(op(i,a),integer) or type(op(i,a),fraction) then
    f:=op(i,a):
  RETURN([f,a-f]):
 fi:
od:
[0,a]:
end:

#Bkna(k,n,a): The pair of rational numbers [c(n),d(n)] such that Akna(k,n,a)=c(n)+d(n)*Akna(k,0,a) of FAIL.
#It also returns the ratio that approximates Akna(k,0,a), the error and the empirical delta
#Try:
#Bkna(3,3,1);
Bkna:=proc(k,n,a) local mu,X,gu,gu1,gu2,ku:
option remember:
gu:=AknaF(k,n,a,X):

mu:=Akna(k,0,a):

gu1:=coeff(gu,X,0):
gu2:=coeff(gu,X,1):

ku:=deltE(-gu1/gu2,mu):

if ku<>FAIL then
 ku:=evalf(ku,10):
fi:

if gu2<>0 then
[gu1,gu2,-gu1/gu2, evalf(-gu1/gu2-mu), ku]:
else
FAIL:
fi:

end:


#EvalCF(L): evaluates the continued fraction
EvalCF:=proc(L) local i:
if not type(L,list) and {seq(type(L[i],integer),i=1..nops(L))}<>{true} then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if nops(L)=1 then
 L[1]:
else
L[1]+1/EvalCF([op(2..nops(L),L)]):
fi:

end:


#ReRa(a,K): inputs a real number approximates it with a rational number. K is a big number indicationg where
#to truncate the continued fraction. Try: ReRa(101/89+10^(-100)); It also returns the error
ReRa:=proc(a,K) local gu,i:
gu:=convert(a,confrac):

for i from 1 to nops(gu) while gu[i]<=gu[1]+10*K do od:


gu:=[op(1..i-1,gu)]:

gu:=EvalCF(gu):
[gu,evalf(gu-a)]:

end:





#Aka1Old(k,a,X): Expressing Akna(k,1,a) in terms of X=Akna(k,0,a); Try:
#Aka1Old(3,4,X);
Aka1Old:=proc(k,a,X) local mu0,mu1,gu,lu:
option remember:

if not (type(k,integer) and  k>=1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if not  type(a,numeric) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


mu0:=expand(Akna(k,0,a)):
mu1:=expand(Akna(k,1,a)):

gu:=RP(mu1):

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


lu:=coeff(evalf(gu[2]/mu0),I,0):

lu:=ReRa(lu,10^10):

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

if abs(lu[2])>10^(-20) then
 RETURN(FAIL):
fi:

if abs(evalf(gu[1]+ lu[1]*mu0-mu1))>10^(-trunc(Digits/2)) then
 FAIL:
else
gu[1]+ lu[1]*X:
fi:


end:



#IsRF(R,x): Is the expression R a rational function of x? Try
#IsRF(x/(x+1),x); IsRat(cos(x),x);
IsRF:=proc(R,x) local d1,d2,i:

d1:=degree(numer(R),x):

if d1=FAIL then
 RETURN(false):
fi:

d2:=degree(denom(R),x):

if d2=FAIL then
 RETURN(false):
fi:

if normal(add(coeff(numer(R),x,i)*x^i,i=0..d1)/add(coeff(denom(R),x,i)*x^i,i=0..d2)-R)<>0 then
RETURN(false):
fi:
true:
end:

#Aka1(k,a,X): Expressing Akna(k,1,a) in terms of X=Akna(k,0,a); Try:
#Aka1(3,4,X);
Aka1:=proc(k,a,X) local mu0,mu1,gu,lu,i,ku:
option remember:


if not (type(k,integer) and  k>=1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if not  type(a,numeric) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if k=1 then
 RETURN(Aka1Old(1,a,X)):
fi:

mu0:=int(1/(1+x^k/a),x):
mu1:=int(x^k*(1-x^k)/(1+x^k/a)^2,x):

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

gu:=0:

for i from 1 to nops(mu1) do
if IsRF(op(i,mu1),x) then
  gu:=gu+op(i,mu1):
fi:
od:


lu:=simplify(normal((mu1-gu)/mu0)):
if not (type(lu,fraction) or type(lu, integer)) then
lprint(lu):
 RETURN(FAIL):
fi:

ku:=lu*X+simplify(subs(x=1,gu)):

mu0:=int(1/(1+x^k/a),x=0..1):
mu1:=int(x^k*(1-x^k)/(1+x^k/a)^2,x=0..1):

if abs(evalf(subs(X=mu0,ku)-mu1))>10^(-9) then
  print(`Something is wrong!`):
  RETURN(FAIL):
fi:

ku:

end:

#AknaF(k,n1,a,X): Fast version of Akna(k,n1,a) using the Almkvist-Zeilberger, where X denotes Akna(k,0,a), i.e.
#int(1/(1+x^k/a),x=0..1);
#AknaF(k,n1,a,X): Fast version of Akna(k,n1,a) using the Almkvist-Zeilberger algoirthm. Try:
#AknaF(3,10,3,X);
#Try:
#AknaF(k,n1,a,X);
AknaF:=proc(k,n1,a,X) local x,ope,n,N,i:
option remember:

if not (type(k,integer) and type(n1,integer) and n1>=0 and k>=1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if not  type(a,numeric) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if Aka1(k,a,X)=FAIL then
  RETURN(FAIL):
fi:

if n1=0 then
 RETURN(X):
fi:
if n1=1 then
 RETURN(Aka1(k,a,X)):
fi:

ope:=AZd1(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x,n,N):


add(subs(n=n1,coeff(ope,N,-i))*AknaF(k,n1-i,a,X),i=1..2):

end:



#ChA(k,a,N0): Checks procedure AknaF(k,n,a,X) for n from 0 to N0. Try:
#ChA(2,2,20);
ChA:=proc(k,a,N0) local n1,mu0,gu1,gu2,X,i1,ku:
mu0:=Akna(k,0,a):
gu1:=[seq(Akna(k,n1,a),n1=0..N0)]:
gu2:=[seq(AknaF(k,n1,a,X),n1=0..N0)]:

if subs(X=mu0,gu2)=gu1 then
 print(`Definitely true`):
 RETURN(true):
else
 print(`The difference in floating point is`):
 ku:=subs(X=evalf(mu0),gu2)-evalf(gu1):
 print(  max(seq(abs(ku[i1]),i1=1..nops(ku)))):
 RETURN(FAIL):
fi:

end:



#c1(a): The constant c1(k,a) such that |Akna(k,n,a)|<=C/c1(k,a)^n. Try:
#c1(a); c1(7,5);
c1:=proc(a) local y,f,ka:

f:=y*(1-y)/(1+y/a):

ka:=numer(normal(diff(f,y)));

ka:=solve(ka,y):

if  0<evalf(ka[1]) and evalf(ka[1])<1 then
ka:=ka[1]:
else
ka:=ka[2]:
fi:

ka:=simplify(1/subs(y=ka,f));
end:


#c2(k,a): The constant c2(a,k) such that |Bkna(k,n,a)[1]|<=C*c2(a,k)^n and |Bkna(k,n,a)[2]|<=C*c2(a,k)^n
#a is a positive integer, k is a pos. integer. Try
#c2(5,3);
c2:=proc(k,a) local ope,n,N,x:
ope:=AZd(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x,n,N)[1]:

ope:=lcoeff(ope,n):
max(solve(ope,N)):
end:


#Seqs(k,a,N): The two sequence B(k,n,a)[1] and B(k,n,a)[2] for n from 1 to N. Try:
#Seqs(3,2,20);
Seqs:=proc(k,a,N) local gu,n1,X:
option remember:


if not  type(a,numeric) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


if Aka1(k,a,X)=FAIL then
  RETURN(FAIL):
fi:

gu:=[seq(AknaF(k,n1,a,X),n1=1..N)]:
[[seq(-coeff(gu[n1],X,0),n1=1..N)],[seq(coeff(gu[n1],X,1),n1=1..N)]]:
end:

#COV(L,K1,K2): inputs a list of rational numbers L, a positive integer K1, and a rational number K2, outputs the sequence L[n]*dn(K1*n)*K2^n1
#Try:
#COV([Seqs(3,2,20)[1],3,1/8);
COV:=proc(L,K1,K2) local n1:
[seq(L[n1]*dn(K1*n1)*K2^n1,n1=1..nops(L))]:
end:


#deltP(k,a,K1,K2): let alpha:=int(1/(1+x^k/a),x=0..1). Write
#int(x^n*(1-x)^n/(1+(x/a)^k)^(n+1),x=0..1)=A(n)-B(n)*alpha
#where A(n) and B(n) are specific rational numbers.
#Assume that you proved that A'(n):=A(n)*lcm(1..K1*n)/K2^n,  B'(n):=B(n)*lcm(1..K1*n)/K2^n, 
#are integers, it outputs the number delta such that
#|alpha-A'(n)/B'(n)|<=CONSTANT/B'(n)^(1+delta). If delta is >0, then this proves the irrationality of
#alpha, and the irrationality measure is 1+1/delta. Try:
#deltP(3,2,2,3);
deltP:=proc(k,a,K1,K2):
(log(c2(k,a))+log(c1(a))) /(log(c2(k,a))+K1-log(K2))
-1:
end:




#Deco11(a,n,r): Given a rational number a, find the smallest K1 such that the denominator of a*dn(n*K1) has <=r factors
Deco11:=proc(a,n,r) local K1:

for K1 from 0 to 20 do
 if  nops(ifactors(denom(a*dn(n*K1)))[2])<=r then
   RETURN([K1,ifactor(denom(a*dn(n*K1)))]    )
 fi:
od:
FAIL:
end:







##############end generalized arctan(1/a)



#GrabHP(n): inputs an integer n and outputs [m,r] such that n/m^r is an integer and r is as large as possible. Try:
#GrabHP(10^20*101);
GrabHP:=proc(n) local gu,r,mu,i:
gu:=ifactors(n)[2]:

r:=max(seq(gu[i][2],i=1..nops(gu))):

mu:=1:

for i from 1 to nops(gu) do
 if gu[i][2]=r then
    mu:=mu*gu[i][1]:
 fi:
od:

[mu,r]:

end:

#GuessDiv(L): Given a sequence of integers L, guesses an integer m, and a rational number such that
#L[n]/m^(trunc(n*alpha)) are integers. Try:
#GuessDiv([seq(5^(3*i)*13,i=1..10)]);

GuessDiv:=proc(L) local i,lu,m:


lu:=[seq(GrabHP(L[i]),i=1..nops(L))]:

if nops({seq(lu[i][1],i=trunc(nops(lu)/2)..nops(lu))})<>1 then
 RETURN(FAIL):
fi:

m:=lu[-1][1]:
m,[seq(i/lu[i][2],i=trunc(nops(lu)/2)..nops(lu))]:

end:


#Rnus(N): The first N positive rational numbers less than 1, Try:
#Rnus(10);
Rnus:=proc(N) local gu,n,a,b:
gu:=[]:

for n from 3 while nops(gu)<=N do
for a from trunc((n-1)/2) to n-1 do
b:=n-a:
if b/a<1 and not member(b/a,{op(gu)}) then
 gu:=[op(gu),b/a]:
fi:
od:
od:

sort([op(1..N,gu)]):

end:








#CheckDH(k,a,N): Checks the divisibility hypothesis for Seqs(k,a,N). Try:
#CheckDH(7,11,100);
CheckDH:=proc(k,a,N) local gu,i,gu1,gu2:
gu:=Seqs(k,a,N);

gu1:=[seq(gu[1][i]*dn(k*i)*k^trunc((k+1)/k*i)/a^i,i=1..N)]:
gu2:=[seq(gu[2][i]*dn(k*i)*k^trunc((k+1)/k*i)/a^i,i=1..N)]:

if evalb(denom(gu1)=[1$nops(gu1)]) and evalb(denom(gu2)=[1$nops(gu2)]) then
 RETURN(true):
else
print([denom(gu1[-1]),denom(gu1[-2])]):
RETURN(false):
fi:


end:

#Tikva(k,a,N1,N2): Given k and a and   positive integers N1 and N2  outputs the
#empirical delta for N1 to N2 for the sequence of rational numbers approximating int(1/(1+x^k/a),x=0..1);
#Try Tikva(1,1,30,40);
Tikva:=proc(k,a,N1,N2) local gu,mu0,i,lu:

if not (type(k,integer) and k>=1 and (type(a,fraction) or type (a,integer)) and type(N1,integer) and N1>=20 and type(N2,integer) and N2>N1) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

 gu:=Seqs(k,a,N2):

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

 mu0:=Akna(k,0,a):


 if abs(evalf(gu[1][N2]/gu[2][N2]-mu0))>10^(-10) then
  print(`Something is wrong`):
  RETURN(FAIL):
 fi:

 if abs(evalf(gu[1][N2]/gu[2][N2]-mu0))<10^(-Digits+200) then
   print(`N is too big, make it smaller`):
   RETURN(FAIL):
 fi:

lu:=[seq(evalf(deltE(gu[1][i]/gu[2][i],mu0)), i=N1..N2)]:

lu:=evalf(lu,10);

min(op(lu)):

end:


#ConjK1(L): inputs a sequence of rational numbers, L, and a limit G, outputs a rational number with numerator
#and denominator <=G such that L[i]*K^i are all integers. Try:
#ConjK1([seq((2/3)^i,i=1..20)],30);

ConjK1:=proc(L,G) local i,K1,K2,L1,K:

if denom(L)=[1$nops(L)] then
 K1:=1:
else

for K1 from 2 to G while denom([seq(L[i]*K1^i,i=1..nops(L))])<>[1$nops(L)] do od:

if K1=G+1 then
 RETURN(FAIL):
fi:

fi:

L1:=[seq(L[i]*K1^i,i=1..nops(L))]:

for K2 from G by -1 to 1  while  {seq(type(L1[i]/K2^i,integer),i=1..nops(L1))}<>{true} do od:

K:=K2/K1:

if not {seq( type(L[i]/K^i,integer), i=1..nops(L)) }<>{integer} then
 RETURN(FAIL):
fi:

K:

end:

#ConjK(k,a,N,G): Given a positive integer k and a rational number a, and a fairly large positive integer N,
#outputs the rational number K with numerators and denominators <=G 
#such that the the sequences of numerators and denominators in the diophantine approximations to int(1/(1+x^k.a),x=0..1) 
#are such that if you multiply the n-th entry by lcm(1..k*n) and divide by K^n you always get an integer, up to the N-th term.
#Try:
#ConjK(2,11,40,30);
ConjK:=proc(k,a,N,G)  local gu,mu1,mu2,i,K1,K2,K1t,K1b,lu:
gu:=Seqs(k,a,N):
if gu=FAIL then
 RETURN(FAIL):
fi:

mu1:=[seq(gu[1][i]*dn(k*i),i=1..N)]:
mu2:=[seq(gu[2][i]*dn(k*i),i=1..N)]:

K1:=ConjK1(mu1,G):

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

K2:=ConjK1(mu2,G):

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

if K1<>K2 then
 print(K1,K2):
 RETURN(FAIL):
fi:




K1t:=numer(K1):
K1b:=denom(K1):

lu:=[[K1b,1],[K1t,1]]:

if K1b=1 then
 RETURN([K1,lu]):
fi:

if denom([seq(gu[1][i]*dn(k*i)*K1b^i/K1t^i,i=1..N)])<>[1$N] or denom([seq(gu[2][i]*dn(k*i)*K1b^i/K1t^i,i=1..N)])<>[1$N] then
 RETURN(FAIL):
fi:

if denom([seq(gu[1][i]*dn(k*i)*K1b^(trunc(i/2))/K1t^i,i=1..N)])=[1$N] and denom([seq(gu[2][i]*dn(k*i)*K1b^(trunc(i/2))/K1t^i,i=1..N)])=[1$N] then
 K1:=K1t/sqrt(K1b):
lu:=[[K1b,1/2],[K1t,1]]:
 RETURN([K1,lu]):
fi:


if denom([seq(gu[1][i]*dn(k*i)*K1b^(trunc(3*i/4))/K1t^i,i=1..N)])=[1$N] and denom([seq(gu[2][i]*dn(k*i)*K1b^(trunc(3*i/4))/K1t^i,i=1..N)])=[1$N] then
 K1:=K1t/K1b^(3/4):
lu:=[[K1b,3/4],[K1t,1]]:
 RETURN([K1,lu]):
fi:



[K1,lu]:

end:


#Info(k,a,G,N): Inputs a pos. integer k, a rational number a, a fairly large (say 100) pos. integer G, and a fairly large
#pos. integer N (say 200), outputs a list [constant,K,deltP,evalf(deltP), deltE] where constant is the  exact value of
#int(1/(1+x^k/a,x=0..1), the constant whose irrationality we are trying to prove,
#K is the rational number such that
#if A(n) B(n) are the sequences of rational number such that A(n)/B(n) is a diaphontine approximation to
#int(1/(1+x^k/a),x=0..1), we have conjecturally that both A(n)*K^n, and B(n)*K^n are always integers (proved, so far for n<=N)
#where the numerator and denominator of K are <=G. deltP is the smallest empirical delta from N/2 to N. Try:
#Info(2,11,100,200);
Info:=proc(k,a,G,N) local K,a1,b1:
 K:=ConjK(k,a,N,G):

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

K:=K[1]:

a1:=deltP(k,a,k,K):
 b1:=Tikva(k,a,trunc(N/2),N):

[Akna(k,0,a),K,a1,evalf(a1,10),b1]:

end:



#Start HMV paper

#kstar(k,a): inputs a pos. integer k and a rational mumber a= -s/r, outputs the denomiator of (s-r)/k 
#and (s-r)/(k*prod(p, p|k). Try:
#kstar(2,-1/7);
kstar:=proc(k,a) local r,s,i,gu1,gu2,p,a1:
a1:=-1/a:
r:=numer(a1):
s:=denom(a1):


if s<0 then
r:=-r:
s:=-s:
fi:

gu1:=denom((s-r)/k):
gu2:=(s-r)/k:

for i from 1 while ithprime(i)<=k do
 p:=ithprime(i):
 if type(k/p,integer) then
  gu2:=gu2/p:
 fi:
od:

gu2:=denom(gu2):

[gu1,gu2]:
end:


#muh(h): the product of p^(1/(p-1)) over all prime divisors of h
muh:=proc(h) local gu,i,p:
gu:=1:

for i from 1 while ithprime(i) <=h do
 p:=ithprime(i):

 if type(h/p,integer) then
  gu:=gu*p^(1/(p-1)):
 fi:
od:
gu:
end:

#lamh(h): h/phi(h)*sum(1/i,i=1..h, gcd(i,h)=1). Try:
#lamh(20);
lamh:=proc(h) local gu1,gu2,i:
gu1:=0:
gu2:=0:
for i from 1 to h do
 if igcd(i,h)=1 then
   gu1:=gu1+1:
    gu2:=gu2+h/i:
 fi:
od:
gu2/gu1:
end:


#IrrM(k,a): The irrationality measure according to HMV with beta=1 . Try 
#IrrM(2,11);
IrrM:=proc(k,a) local a1,r,s,gu,gu3:


a1:=-1/a:
r:=numer(a1):
s:=denom(a1):

if s<0 then
r:=-r:
s:=-s:
fi:

print(`r is`, r):
print(`s is`, s):

gu:=kstar(k,a):

gu3:=log(min(gu[1]*muh(k),gu[2])):

1-
(2*log(sqrt(s)+sqrt(s-r))+lamh(k)+gu3)/
(2*log(abs(sqrt(s)-sqrt(s-r)))+lamh(k)+gu3):
end:




#End HMV paper


##start Story procedures


#Paper(k,a,G,N): A paper about the irrationality of int(1/(1+x^k/a),x=0..1); G and N are parameters for estimating and guessing. Try:
#Paper(2,11,100,200):
Paper:=proc(k,a,G,N) local c,gu,x,E,n,A,B,ope,SN,X,lu,vu,d,y,delta, C, i:

gu:=Info(k,a,G,N):

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

if gu[4]<0 then
 print(`Can't prove irrationality for`, Int(1/(1+x^k/a),x=0..1)):
 RETURN(FAIL):
fi:

print(`Proof of Irrationaility, and an Irrationality Measure of`,   Int(1/(1+x^k/a),x=0..1)):
print(``):
print(`Proposition: Let c be the constant`,  Int(1/(1+x^k/a),x=0..1), `that happens to be equal to`, gu[1], ` alias`, evalf(gu[1],20)):
print(``):
print(`Then c is irrational, and has irrationality measure at most`, simplify(normal(1+1/gu[3])), `that equals`, evalf(1+1/gu[3],20)):
print(``):
print(`Proof: Consider `):
print(``):
print(E(n)=Int(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x=0..1)):
print(``):
print(`This can be written as `):
print(``):
print(E(n)=B(n)*c-A(n) ):
print(``):
print(`For some sequences`, A(n), B(n), `Of RATIONAL NUMBERS `):
print(``):
print(`By looking at the maximum of`, y*(1-y)/(1+y/a), ` in 0<=y<=1, (note that`, y=x^k, ` ), it is readily seen that`):
print(``):
print(E(n)<=C/c1(a)^n):
print(``):
ope:=AZd1(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x,n,SN):

print(``):
print(`It follows from the  Almkvist-Zeilberger algorithm that  E(n) satisfies the recurrence `):
print(``):
print(E(n)=add(coeff(ope,SN,-i)*E(n-i),i=1..-ldegree(ope,SN))):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(E(0)=c , E(1)=AknaF(k,1,a,c)):


print(`and in Maple format `):
print(``):
lprint(E(n)=add(coeff(ope,SN,-i)*E(n-i),i=1..-ldegree(ope,SN))):
print(``):
lprint(`Subject to the initial conditions`):
print(``):
lprint(E(0)=c , E(1)=AknaF(k,1,a,c)):


print(``):
print(`It follows that the sequence of rational numbers`, A(n), B(n), `satisfy the same recurrence `):
print(``):
print(B(n)=add(coeff(ope,SN,-i)*B(n-i),i=1..-ldegree(ope,SN))):
print(``):
print(`subject to the intial conditions`):
print(``):
print(B(0)=1, B(1)=coeff(AknaF(k,1,a,X),X,1)):
print(``):


print(`and in Maple format `):

print(``):
print(`It follows that the sequence of rational numbers`, A(n), B(n), `satisfy the same recurrence `):
print(``):
lprint(B(n)=add(coeff(ope,SN,-i)*B(n-i),i=1..-ldegree(ope,SN))):
print(``):
lprint(`subject to the intial conditions`):
print(``):
lprint(B(0)=1, B(1)=coeff(AknaF(k,1,a,X),X,1)):
print(``):


print(` and `):
print(``):
print(A(n)=add(coeff(ope,SN,-i)*A(n-i),i=1..-ldegree(ope,SN))):
print(``):
print(`subject to the intial conditions`):
print(``):
print(A(0)=0, A(1)=coeff(AknaF(k,1,a,X),X,0)):
print(``):
print(`The sequence of rational numbers`, A(n)/B(n), `are good numerical approximations for c`):




print(` In Maple format: `):
print(``):
lprint(A(n)=add(coeff(ope,SN,-i)*A(n-i),i=1..-ldegree(ope,SN))):
print(``):
lprint(`subject to the intial conditions`):
print(``):
lprint(A(0)=0, A(1)=coeff(AknaF(k,1,a,X),X,0)):
print(``):
print(``):

print(`The sequence of rational numbers`, A(n)/B(n), `are good numerical approximations for c`):


print(``):


lu:=Seqs(k,a,N):
if lu=FAIL then
 RETURN(FAIL):
fi:

print(`Just for fun`):

print(` The `, N , `-th term of that sequence is`):
print(``):
print(lu[1][N]/lu[2][N]):
print(``):
print(`and its differene from  c is`):
print(``):
print(evalf(abs(lu[1][N]/lu[2][N]-gu[1]))):
print(``):
print(`The smallest empirical delta from`, N/2, ` to `, N, `is `):
print(``):
print(gu[5]):
print(``):
print(`Alas, the sequences A(n), and B(n) are not integers, but rational mumbers. We need the following divisibility lemma`):
print(``):
print(`that we leave to the reader `):
print(``):
vu:=ConjK(k,a,N,G)[2]:

print(`Let d(n) be the least common multiple of the first n natural numbers, 1...n `):

if vu[1][2]=1 then
 print(`Lemma: ` , A(n)*d(k*n)*vu[1][1]^n/vu[2][1]^n, `and `,  B(n)*d(k*n)*vu[1][1]^n/vu[2][1]^n, ` are always integers `):
else
 print(`Lemma: ` , A(n)*d(k*n)*vu[1][1]^(trunc (n*vu[1][2]) )  /vu[2][1]^n, `and `,  B(n)*d(k*n)*vu[1][1]^(trunc (n*vu[1][2]))/vu[2][1]^n, ` are always integers `):
fi:

print(`Let's  call these new integer sequences`, A1(n), B1(n) ):
print(``):
print(`Using the above recurrence, by the Poincare lemma, the rate of growth of A(n), B(n) are (essentially)`):
print(``):
print(C*c2(k,a)^n):
print(``):
print(`Dividing `, E(n)=B(n)*c-A(n), ` by `, B(n), ` we get `):
print(``):
print(abs(c-A(n)/B(n)) <= C/(c1(a)*c2(a,k))^n):
print(``):	       
print(`Hence `):
print(``):
print(abs(c-A1(n)/B1(n)) <= C/(c1(a)*c2(a,k))^n):
print(``):	       
print(`But `, B1(n)=B(n)*d(k*n)*vu[1][1]^(trunc (n*vu[1][2]))/vu[2][1]^n, `hence `):
print(``):
vu:=ConjK(k,a,N,G)[1]:
print(``):
print(B1(n), `is of the order `, (exp(k)*c2(a,k)*vu)^n):
print(``):
print(`Hence `):
print(``):
print(abs(c-A1(n)/B1(n)) ):
print(` is <=  `, C/B1(n)^(1+delta)):
print(``):
print(`where delta equals`, gu[3]):
print(``):
print(`That in floating-point is`, gu[4]):
print(``):
print(`It follows that an irrationality measure for c is`):
print(``):
print(simplify(normal(1+1/gu[3]))):
print(``):
print(`that equals, approximately`):
print(``):
print(evalf(1+1/gu[4],10)):
print(``):
print(`------------------------------------------`):

end:






#Prop1(k,a,G,N,co): A proposition numbered co, about the irrationality of int(1/(1+x^k/a),x=0..1); G and N are parameters for estimating and guessing. Try:
#Prop1(2,11,100,200,10):
Prop1:=proc(k,a,G,N,co) local c,gu,x,E,n,A,B,ope,SN,X,lu,vu,d,delta, C, i, y:


gu:=Info(k,a,G,N):

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

if gu[4]<0 then
 print(`Can't prove irrationality for`, Int(1/(1+x^k/a),x=0..1)):
 RETURN(FAIL):
fi:


print(`Proposition Number`, co, `: Let c be the constant`,  Int(1/(1+x^k/a),x=0..1), `that happens to be equal to`, gu[1], ` alias`, evalf(gu[1],20)):
print(``):
print(`Then c is irrational, and has irrationality measure at most`, simplify(normal(1+1/gu[3])), `that equals`, evalf(1+1/gu[3],20)):
print(``):
print(`Proof: Consider `):
print(``):
print(E(n)=Int(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x=0..1)):
print(``):
print(`This can be written as `):
print(``):
print(E(n)=B(n)*c-A(n) ):
print(``):
print(`For some sequences`, A(n), B(n), `Of RATIONAL NUMBERS `):
print(``):
print(`By looking at the maximum of`, y*(1-y)/(1+y/a), ` in 0<=y<=1, (note that`, y=x^k, `), it is readily seen that`):
print(``):
print(E(n)<=C/c1(a)^n):
print(``):
ope:=AZd1(x^(k*n)*(1-x^k)^n/(1+x^k/a)^(n+1),x,n,SN):

print(``):
print(`It follows from the  Almkvist-Zeilberger algorithm that  E(n) satisfies the recurrence `):
print(``):
print(E(n)=add(coeff(ope,SN,-i)*E(n-i),i=1..-ldegree(ope,SN))):
print(``):
print(`Subject to the initial conditions`):
print(``):
print(E(0)=c , E(1)=AknaF(k,1,a,c)):


print(`and in Maple format `):
print(``):
lprint(E(n)=add(coeff(ope,SN,-i)*E(n-i),i=1..-ldegree(ope,SN))):
print(``):
lprint(`Subject to the initial conditions`):
print(``):
lprint(E(0)=c , E(1)=AknaF(k,1,a,c)):


print(``):
print(`It follows that the sequence of rational numbers`, A(n), B(n), `satisfy the same recurrence `):
print(``):
print(B(n)=add(coeff(ope,SN,-i)*B(n-i),i=1..-ldegree(ope,SN))):
print(``):
print(`subject to the intial conditions`):
print(``):
print(B(0)=1, B(1)=coeff(AknaF(k,1,a,X),X,1)):
print(``):


print(`and in Maple format `):

print(``):
print(`It follows that the sequence of rational numbers`, A(n), B(n), `satisfy the same recurrence `):
print(``):
lprint(B(n)=add(coeff(ope,SN,-i)*B(n-i),i=1..-ldegree(ope,SN))):
print(``):
lprint(`subject to the intial conditions`):
print(``):
lprint(B(0)=1, B(1)=coeff(AknaF(k,1,a,X),X,1)):
print(``):


print(` and `):
print(``):
print(A(n)=add(coeff(ope,SN,-i)*A(n-i),i=1..-ldegree(ope,SN))):
print(``):
print(`subject to the intial conditions`):
print(``):
print(A(0)=0, A(1)=coeff(AknaF(k,1,a,X),X,0)):
print(``):
print(`The sequence of rational numbers`, A(n)/B(n), `are good numerical approximations for c`):




print(` In Maple format: `):
print(``):
lprint(A(n)=add(coeff(ope,SN,-i)*A(n-i),i=1..-ldegree(ope,SN))):
print(``):
lprint(`subject to the intial conditions`):
print(``):
lprint(A(0)=0, A(1)=coeff(AknaF(k,1,a,X),X,0)):
print(``):
print(``):

print(`The sequence of rational numbers`, A(n)/B(n), `are good numerical approximations for c`):


print(``):


print(`Just for fun`):
lu:=Seqs(k,a,N):
print(` The `, N, ` -th term of that sequence is`):
print(``):
print(lu[1][N]/lu[2][N]):
print(``):
print(`and its differene from  c is`):
print(``):
print(evalf(abs(lu[1][N]/lu[2][N]-gu[1]))):
print(``):
print(`The smallest empirical delta from`, N/2, ` to `, N, `is `):
print(``):
print(gu[5]):
print(``):
print(`Alas, the sequences A(n), and B(n) are not integers, but rational mumbers. We need the following divisibility lemma`):
print(``):
print(`that we leave to the reader `):
print(``):
vu:=ConjK(k,a,N,G)[2]:

print(`Let d(n) be the least common multiple of the first n natural numbers, 1...n `):

if vu[1][2]=1 then
 print(`Lemma: ` , A(n)*d(k*n)*vu[1][1]^n/vu[2][1]^n, `and `,  B(n)*d(k*n)*vu[1][1]^n/vu[2][1]^n, ` are always integers `):
else
 print(`Lemma: ` , A(n)*d(k*n)*vu[1][1]^(trunc (n*vu[1][2]) )  /vu[2][1]^n, `and `,  B(n)*d(k*n)*vu[1][1]^(trunc (n*vu[1][2]))/vu[2][1]^n, ` are always integers `):
fi:

print(`Let's  call these new integer sequences`, A1(n), B1(n) ):
print(``):
print(`Using the above recurrence, by the Poincare lemma, the rate of growth of A(n), B(n) are (essentially)`):
print(``):
print(C*c2(k,a)^n):
print(``):
print(`Dividing `, E(n)=B(n)*c-A(n), ` by `, B(n), ` we get `):
print(``):
print(abs(c-A(n)/B(n)) <= C/(c1(a)*c2(a,k))^n):
print(``):	       
print(`Hence `):
print(``):
print(abs(c-A1(n)/B1(n)) <= C/(c1(a)*c2(a,k))^n):
print(``):	       
print(`But `, B1(n)=B(n)*d(k*n)*vu[1][1]^(trunc (n*vu[1][2]))/vu[2][1]^n, `hence `):
print(``):
vu:=ConjK(k,a,N,G)[1]:
print(``):
print(B1(n), `is of the order `, (exp(k)*c2(a,k)*vu)^n):
print(``):
print(`Hence `):
print(``):
print(abs(c-A1(n)/B1(n)) ):
print(` is <=  `, C/B1(n)^(1+delta)):
print(``):
print(`where delta equals`, gu[3]):
print(``):
print(`That in floating-point is`, gu[4]):
print(``):
print(`It follows that an irrationality measure for c is`):
print(``):
print(simplify(normal(1+1/gu[3]))):
print(``):
print(`that equals, approximately`):
print(``):
print(evalf(1+1/gu[4],10)):
print(``):
print(`------------------------------------------`):

end:







#Book(k, Maxa, G,N): Inputs a positive integer k, a positive integer Maxa, and positive integers G and N (parameters for guessing)
#outputs a book about irrationaility, and irrationality measures of int(1/(1+x^k/a),x=0..1) for a from 1 to Maxa
#Try:
#Book(1,20,200,200);
Book:=proc(k, Maxa, G,N) local a,co,gu,x,t0:
t0:=time():
print(``):
print(`On the Irrationality of`, Int(1/(1+x^k/a),x=0..1), `for a from 1 to `,  Maxa):
print(``):
print(` By Shalosh B. Ekhad`):
print(``):
print(`In this computer-generated book, accompanying the article by Doron Zeilberger and Wadim Zudilin`):
print(` "Towards Automatic Discovery of Irrationailty Proofs and Irrationality Measures."`):
print(``):
print( ` I will automatically prove irrationality, and establish irrationality measures, for the constant`):
print(``):
print(Int(1/(1+x^k/a),x=0..1)):
print(``):
print(`for a from 1 to`, Maxa):
print(``):
print(`I will also state the cases where we were unable to do it.`):
print(``):

co:=0:

for a from 1 to Maxa do
 gu:=Info(k,a,G,N):

if gu[5]<0 then
 print(`We are unable, with this method to prove the irrationality of`, Int(1/(1+x^k/a),x=0..1)):
elif gu[4]<0 then
 print(`We are unable, with this method to prove the irrationality of`, Int(1/(1+x^k/a),x=0..1), `but it seems that a sharperning of`):
 print(` the divisibility lemma might lead to a proof`):
print(`since the smallest empirical delta from`, N/2, `to `, N, ` is positive:`, gu[5]):
else
co:=co+1:
print(``):
print(`-----------------------------------------------`):
Prop1(k,a,G,N,co):
print(``):
fi:

od:

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