#####################################################################
##CatC.txt: Save this file as  CatC.txt                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  CatC.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.30, 2019`):
print(` This is CatC.txt `):
print(`It is one of the Maple packages that accompany the article `):
print(` Automatic Discovery of Irrationailty Proofs and Irrationality Measures`):
print(`by Doron Zeilberger and Wadim Zudiln`):

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



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

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

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

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

with(combinat):



ezra1:=proc()

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

else
ezra(args):
fi:

end:

ezra:=proc()

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


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



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



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



elif nops([args])=1 and op(1,[args])=MakeNice then
print(`MakeNice(ope,n,N): inputs a recurrence operator ope in n, N, outputs Ope1(n,1/N) such that`):
print(`1-Ope1(n,N) is equivalent fo Ope. Try:`):
print(`MakeNice(N^2-n*N-1,n,N);`):

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

elif nops([args])=1 and op(1,[args])=Paper then
print(`Paper():  A computer generated paper about the joint approximation of 1, Catalan, Pi, and log(2). Try:`):
print(`Paper(): `):

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



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


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

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


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


elif nops([args])=1 and op(1,[args])=Z2gF then
print(`Z2gF(n): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1+x^2*y^2)^(n+1),x=0..1),y=0..1);`):
print(`Done FAST, using the fourth-order linear recurrence satisfied by it, found by the multi-Almkvist-Zeilberger`):
print(`algorithm. Try`):
print(`[seq(Z2gF(i),i=0..10)];`):


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

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

Digits:=7000:

#########From EKHAD

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

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

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


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

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

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

0:
end:

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


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



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

Digits:=3000:






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



#Z2gF(n): int(int(x^n*(1-x)^n*y^n*(1-y)^n/(1+x^2*y^2)^(n+1),x=0..1),y=0..1) The fast way, using the recurrence
Z2gF:=proc(n)  option remember:
if n=0 then
Catalan:
elif n=1 then
 -1+1/2*Catalan+1/8*Pi+1/4*ln(2)
elif n=2 then
 -1+1/2*ln(2)+1/16*Pi+1/2*Catalan
elif n=3 then
 -185/144+19/24*ln(2)+7/8*Catalan-1/48*Pi
else

(129*n^4-663*n^3+1100*n^2-710*n+162)/(43*n^2-178*n+171)/n^2*Z2gF(n-1)
-1/2*(215*n^4-1320*n^3+2848*n^2-2537*n+756)/(43*n^2-178*n+171)/n^2*Z2gF(n-2)
+1/4*(387*n^4-2763*n^3+6896*n^2-\6838*n+2052)/(43*n^2-178*n+171)/n^2*Z2gF(n-3)
-1/8*(43*n^2-92*n+36)*(n-3)^2/(43*n^2-178*n+171)/n^2*Z2gF(n-4)
fi:

end:


#Z2gFS(n,X1,X2,X3): Like Z2gF(n) but in terms of X1, X2,X3,  standing for  Catalan, Pi, and log(2), respecitvely. Try:
#Z2gFS(10,3,X1,X2,X3);
Z2gFS:=proc(n,X1,X2,X3)  option remember:

if n=0 then
X1:
elif n=1 then
 -1+1/2*X1+1/8*X2+1/4*X3
elif n=2 then
 -1+1/2*X3+1/16*X2+1/2*X1
elif n=3 then
 -185/144+19/24*X3+7/8*X1-1/48*X2
else
(129*n^4-663*n^3+1100*n^2-710*n+162)/(43*n^2-178*n+171)/n^2*Z2gFS(n-1,X1,X2,X3)
-1/2*(215*n^4-1320*n^3+2848*n^2-2537*n+756)/(43*n^2-178*n+171)/n^2*Z2gFS(n-2,X1,X2,X3)
+1/4*(387*n^4-2763*n^3+6896*n^2-\6838*n+2052)/(43*n^2-178*n+171)/n^2*Z2gFS(n-3,X1,X2,X3)
-1/8*(43*n^2-92*n+36)*(n-3)^2/(43*n^2-178*n+171)/n^2*Z2gFS(n-4,X1,X2,X3)
fi:

end:








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

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

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


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

end:




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

gu:=min(seq(deltEg(numer(Z2gFS(n,X1,X2,X3)),X1,X2,X3,Catalan,Pi,log(2)),n=N..2*N)):
evalf(gu,20):
end:



#Ope(n,N): The recurrence operator annihilating  int(x^n*(1-x)^n*y^n*(1-y)^n/(1+x^2*y^2)^(n+1). Try:
#Ope(n,N);
Ope:=proc(n,N)
1/8*(43*n^2+252*n+356)*(1+n)^2/(43*n^2+166*n+147)/(n+4)^2-1/4*(387*n^4+3429*n^3
+10892*n^2+14778*n+7276)/(43*n^2+166*n+147)/(n+4)^2*N+1/2*(215*n^4+2120*n^3+
7648*n^2+11927*n+6736)/(43*n^2+166*n+147)/(n+4)^2*N^2-(129*n^4+1401*n^3+5528*n^
2+9290*n+5514)/(43*n^2+166*n+147)/(n+4)^2*N^3+N^4:
end:


#MakeNice(ope,n,N): inputs a recurrence operator ope in n, N, outputs Ope1(n,1/N) such that
#1-Ope1(n,N) is equivalent fo Ope. Try:
#MakeNice(N^2-n*N-1,n,N);
MakeNice:=proc(ope,n,N) local k,ope1:
k:=degree(ope,N):
ope1:=expand(subs(n=n-k,ope)/N^k):
ope1:=1-ope1:

add(factor(coeff(ope1,N,-i))*N^(-i),i=1..k):
end:

Ope1:=proc(N) local n:
factor(lcoeff(numer(Ope(n,N)),n)):
end:


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







#ZZg(): The smallest and largest roots of the indicial equation
ZZg:=proc() local gu,N,ka,ga,i,si:
gu:=evalf([solve(Ope1(N),N)]):

ka:=gu[1]:
si:=evalf(gu[1]):
for i from 2 to nops(gu) do
 if abs(gu[i])<ka then
   ka:=gu[i]:
   si:=abs(ka):
 fi:
od:


ga:=gu[1]:
si:=abs(gu[1]):
for i from 2 to nops(gu) do
 if abs(gu[i])>ga then
   ga:=gu[i]:
   si:=abs(ga):
 fi:
od:
[ka,ga]:

end:


#deltT(): The theoretical delta such that |A(n)+B(n)*Catalan+C(n)*Pi+ D(n)*log(2))|<=C/max(A(n),B(n),C(n),D(n))^delta. Try
#deltT();
deltT:=proc() local gu,N,ka,ga,i,si:
gu:=evalf([solve(Ope1(N),N)]):

ka:=gu[1]:
si:=evalf(gu[1]):
for i from 2 to nops(gu) do
 if abs(gu[i])<ka then
   ka:=gu[i]:
   si:=abs(ka):
 fi:
od:


ga:=gu[1]:
si:=abs(gu[1]):
for i from 2 to nops(gu) do
 if abs(gu[i])>ga then
   ga:=gu[i]:
   si:=abs(ga):
 fi:
od:

evalf(-(log(ka)+2+log(4))/(log(ga)+2+log(4)),20):
end:







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

lu:=evalf(ZZg(),20):

print(`Integer Linear Combinations of Catalan,Pi, and log(2)  that are close to an INTEGER (but not close enough) [A FAILED ATTEMPT]`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: There exist EXPLICIT integer sequences A(n), B(n), C(n), D(n) `):
print(`such that `);
print(``):
print(`abs(A(n)+B(n)*Catalan+C(n)*Pi+D(n)*log(2))<= CONSTANT/max(A(n),B(n),C(n),D(n))^delta, where`):
print(`delta equals `):
print(``):
print(gu):
print(``):
print(`Comment 1: Alas, since`, gu, `is negative, this is disappointing. Nevertheless`, gu, `is not that negative and this raises the hope that`):
print(`with other integrals it may work out to have a positive delta, like we did for Dilog((a-1)/a) and Log((a-1)/a) for a integers. `):
print(``):

vu:=Z2gDelE(200):

vu:=evalf(vu,20):

print(`Comment 2: Just for comparison, the smallest empirical delta between 200 and 400 is`, vu):
print(``):
print(`--------------------------------------------------------------------------------`):
print(``):

print(`Proof of the Theorem`):
print(``):
print(`Consider an analog of Beukers's famous integral that established the improved Apery proof of the irrationality of Zeta(2)`):
print(``):
print(`Let's call it`, E(n) ):
print(``):
print(E(n)=Int(Int(x^n*(1-x^n)*y^n*(1-y)^n/(1+x^2*y^2)^(n+1),x=0..1),y=0..1)):
print(``):
print(`It is readily seen that E(n) is a linear combination with RATIONAL number coefficients of 1, Catalan, Pi, and log(2) `):
print(``):
print(`E(n)=A1(n)+B1(n)*Catalan+ C1(n)*Pi+ D1(n)*log(2) `):
eprint(``):
print(`This introduces four sequences of RATIONAL numbers, A1(n), B1(n), C1(n), D1(n) `):
print(``):
print(`Lemma 1: The absolute value of E(n)  is asymptotic (up to a constant) to`):
print(``):
print(lu[1]^n):
print(``):
print(`Proof: Maximize the function`, x*(1-x)*y*(1-y)/(1+x^2*y^2), `in 0<=x,y<=1 using  two-variable calculus`):
print(``):
print(`Indeed taking partial derivatives with respect to x and y, and setting it equal to zero gives that the maximum is attained at the point`):
f:= x*(1-x)*y*(1-y)/(1+x^2*y^2):
ku:=solve({numer(diff(f,x)),numer(diff(f,y))},{x,y})[7]:
print(ku):
print(``):
print(`Now plug it into`,  x*(1-x)*y*(1-y)/(1+x^2*y^2) ):
print(``):
print(`Lemma 2: The four sequences of rational numbers A1(n), B1(n), C1(n), D1(n) are all asymptotic, up to a constant to`):
print(``):
print(lu[2]^n):
print(``):
print(`Proof: Using the amazing Multivariable extension of the Almkvist-Zeilberger algorithm, due to Moa Apagodu and Doron Zeilberger`):
print(`it was discovered and proved that E(n) satisfies the following FOURTH-ORDER linear recurrence  equation with polynomial coefficients`):
vu:=Ope(n,N):
print(``):
print(add(coeff(vu,N,i)*E(n+i),i=0..degree(vu,N))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(add(coeff(vu,N,i)*E(n+i),i=0..degree(vu,N))=0):
print(``):
print(`It follows that the four sequences of rational numbers, A1(n), B1(n), C1(n), D1(n) also satisy this recurrence, so let denote them`):
print(`all by X(n)`):
print(``):
print(add(coeff(vu,N,i)*X(n+i),i=0..degree(vu,N))=0):
print(``):

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