###################################################################### ##SALIKHOVpi.txt: Save this file as SALIKHOVpi.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read SALIKHOVpi.txt # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #DoronZeil at gmail dot com # ###################################################################### #Created: Nov.-Dec. 2019 read `SalikhovDataFile.txt`: Digits:=20000: print(`Created: Nov.-Dec. 2019`): print(` This is SALIKHOVpi.txt `): print(`It is the Maple package that accompanies the article `): print(` The Irrationality Measure of Pi is at most 7.103205334137... `): print(`by Doron Zeilberger and Wadim Zudiln`): print(`available from the arxiv and from `): print(` http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/pimeas.html `): print(``): print(`It is based on V. Kh. Salikhov's beautiful article "On the Measure of Irrationality of the Number Pi"`): print(` Mathematical Notes,2010, vol.88, no. 4, 563-573, but uses an alternative approach `): 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 procedures about the A=2 B=3 case, type ezra23();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): print(`For a list of the procedures about arctan(1/7), type ezraAT7();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): 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): ezraAT7:=proc() if args=NULL then print(` The arctan(1/7) procedures are: AT7Seq, AT7SeqF, PairSeqAT7 `): print(` `): print(``): else ezra(args): fi: end: ezra23:=proc() if args=NULL then print(` The procedures about the A=2, B=3 case are: Info23, Ope23, Paper23, PiSeq23`): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(` The supporting procedures are: AppSeqPi, J, JF, JFSh, JSh, JFShFe, GuessK, GuessK1, GuessKr, MakeDB`): print(` Nu1, OpePiAB, OpeSh, OpeShe, OpeSho, PairSeqPi, PiSeq, PiSeqF, SeqFromRec, `): print(` `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: AppSeqPiAB, AZd, BestAB, EmpDelAB, GuessKrPair, InfoAB, NorPair, Nu1Extra, PairSeqPiAB, PaperAB, PiSeqAB, ShoreshPiAB `): print(` `): elif nops([args])=1 and op(1,[args])=AppSeqAT7 then print(`AppSeqAT7(N): the first N terms of the approximating sequence to arctan(1/7), using the Salikhov sequence. Try:`): print(`AppSeqAT7(20);`): elif nops([args])=1 and op(1,[args])=AppSeqPi then print(`AppSeqPi(N): the first N terms of the Salikhov approximating sequence to Pi.`): print(`For a faster version, use: AppSeqPiAB(3,5,N);`): print(` Try:`): print(`AppSeqPi(20);`): elif nops([args])=1 and op(1,[args])=AppSeqPiAB then print(`AppSeqPiAB(A,B,N): the first N terms of the generalized Salikhov approximating sequence to Pi using SeqPiAB(A,B,N). Try:`): print(`AppSeqPiAB(3,5,100);`): elif nops([args])=1 and op(1,[args])=AT7Seq then print(`AT7Seq(N): the first N terms in the Salikhov odd sequence, featuring arctan(1/7). VERY SLOW, using the definition. For checking purposes only.`): print(`Try:`): print(`AT7Seq(10);`): elif nops([args])=1 and op(1,[args])=AT7SeqF then print(`AT7SeqF(N): the first N terms in the Salikhov odd sequence, featuring arctan(1/7). FAST version of AT7Seq(N),`): print(`using the recurrence gotten by the Almkvist-Zeilberger algorithm.`): print(`Try:`): print(`AT7Seq(10);`): 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])=BestAB then print(`BestAB(K,M): Inputs a positive integer K <=10, and a large positive integer, M`): print(`Explores the empirical delta for the terms between n=M-10 and n=M, in the rational approximations in the generalized Salikhov integral. Using`): print(`Int(1/I*(x-4+2*I)^(A*n1)*(x-4-2*I)^(A*n1)*(x-5)^(A*n1)*(x-6+2*I)^(A*n1)*(x-6-2*I)^(A*n1)/(x^(B*n1+1)*(x-10)^(B*n1+1)),x=4-2*I..4+2*I)`): print(`for all 1<=A,B<=K with gcd(A,B)=1 and (8,13).Try:`): print(`BestAB(10,300);`): elif nargs=1 and args[1]=EmpDelAB then print(`EmpDelAB(A,B,M): the empirical deltas from M-10 to M for the Diophatine approximations implied by`): print(`J(4,2,5,6,10,A,B,2*n0). Try:`): print(`EmpDelAB(1,1,40);`): elif nargs=1 and args[1]=GuessK then print(`GuessK(L,m,r): inputs a sequence of rational numbers, L, a positive integer m, and a positive integer r`): print(` finds a rational number K such that L[i]*K^i*dn(r*i) are all integers. Try:`): print(`GuessK([seq((2/3)^i/dn(5*i),i=1..20)],30,5);`): elif nargs=1 and args[1]=GuessK1 then print(`GuessK1(L,m): Inputs a sequence of integers L, `): print(` finds a constant K of the form K=ithprime(i)^j, with i and j<=m such that L[n]/K^n is an integer for all n. Try:`): print(`GuessK1([seq((i^3+1)*25^(i),i=1..30)],10);`): elif nargs=1 and args[1]=GuessKr then print(`GuessKr(L,m): inputs a sequence of rational numbers, L, a positive integer m,`): print(`finds a positive integer r, and a rational number K [with search space given by m, see GuessK1(L,m)], `): print(`such that L[i]*K^i*dn(r*i) are all integers. Try:`): print(`GuessKr([seq((2/3)^i/dn(5*i),i=1..20)],30);`): elif nargs=1 and args[1]=GuessKrPair then print(`GuessKrPair(P,m): Given a pair,P, of sequences of rational numbers and a positive integer m`): print(`finds a non-negative integer r and a rational number K such that`): print(`BOTH P[1][n]*dn(r*n)*K^(i*n) and P[1][n]*dn(r*n)*K^(i*n) are sequences of INTEGERS. `): print(`It tries out all the primes up to m and all the powers up to m.`): print(`Type: `): print(`GuessKrPair(PairSeqPiAB(3,5,40),20);`): elif nargs=1 and args[1]=InfoAB then print(`InfoAB(A,B,x,n,N,M): inputs positive integers A and B such that gcd(A,B), symbols x,n,N, and a large positive integer M`): print(`outputs a list consisting of `): print(`(i)1/I*Int((x-4+2*I)^(2*A*n)*(x-4-2*I)^(2*A*n)*(x-5)^(2*A*n)*(x-6+2*I)^(2*A*n)*(x-6-2*I)^(2*A*n)/(x^(2*B*n+1)*(x-10)^(2*B*n+1)),x=4-2*I..4+2*I):`): print(`(this is the generalized Salikhov integral)that can be written as C(n)+D(n)*Pi for some rational numbers C(n) and D(n)`): print(`(ii) the Empirical delta for the approximation to Pi from the sequence obtained from the generalized Salikhov integral`): print(`with A, B i.e.`): print(`(iii) the linear recurrence operator in n and the shift operator N annihilating E(n) and hence the coefficients of Pi and 1`): print(`(iv) The leading coefficient in n, in other words the constant-coefficient operator approximinating it.`): print(`the so-called indicial equation`): print(`(v) The pair [largest root, absolute value of smaller roots]`): print(`(vi) The conjectured pair [r,K] such that lcm(1..r*n)*K^n*E(n) has integer coefficients`): print(`In other words C1(n):=lcm(1..r*n)*K^n*C(n) and D1(n):=lcm(1..r*n)*K^n*D(n) are INTEGERS sequences`): print(`(vii) The upper bound for the irrationality measure not taking advantage of the fact that gcd(C1(n),D1(n)) may grow. `): print(`(viii) the smallest log(gcd(A1(n),A2(n))/n for n between M/2 and M`): print(`(ix) assuming that this is a lower bound for the lim sup of log(gcd(C1(n),D1(n))/n as n goes to infinity ,`): print(`the implied upper bound for the irrationality measure of Pi. Try:`): print(`InfoAB(3,5,x,n,N,400);`): elif nargs=1 and args[1]=Info23 then print(`Info23(x,n,N,M): inputs symbols x,n,N, and a large positive integer N, and outputs a list consisting of `): print(`(i)1/I*Int((x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*A*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*A*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x=4-2*I..4+2*I):`): print(`(this is the generalized Salikhov integral) that can be written as C(n)+D(n)*Pi for some rational numbers C(n) and D(n)`): print(`(ii) the Empirical delta for the approximation to Pi from the sequence obtained from the modified Salikhov integral`): print(`with 3 replaced by 2 and 5 replaced by 3, but vaild for all n, not just even n.`): print(`(iii) the linear recurrence operator in n and the shift operator N annihilating E(n) and hence the coefficients of Pi and 1`): print(`(iv) The leading coefficient in n, in other words the constant-coefficient operator approximinating it.`): print(`the so-called indicial equation`): print(`(v) The pair [largest root, absolute value of smaller roots]`): print(`(vi) The conjectured pair [r,K] such that lcm(1..r*n)*K^n*E(n) has integer coefficients`): print(`Let's call the coefficients if 1 and Pi A1(n) and A2(n) respectively`): print(`(vii): The upper bound for the irrationality NOT taking advantage of the fact that gcd(C(n),D(n)) may`): print(`(vii) the smallest log(gcd(A1(n),A2(n))/n for n between M/2 and M`): print(`(viii) assuming that this is a lower bound for the lim sup of log(gcd(A1(n),A2(n))/n for n between M/2 and M`): print(`the implied upper bound for the irrationality measure of Pi. Try:`): print(`Info23(x,n,N,400);`): elif nargs=1 and args[1]=J then print(`J(a1,b1,c1,d1,e1,,A,B,n1): `): print(`1/I*int((x-a1+b1*I)^(A*n)*(x-a1-b1*I)^(A*n)*(x-c1)^(A*n)*(x-d1+b1*I)^(A*n)*(x-d1-b1*I)^(A*n)/(x^(B*n)*(x-e1)^(B*n)),x=a1-b1*I..a1+b1*I):`): print(`The original Salikhov integral is`): print(`J(4,2,5,6,10,3,5,n1);`): elif nargs=1 and args[1]=JF then print(`JF(a1,b1,c1,d1,e1,A,B,n1,x) : the integrand of the generalized Salikhov integral. Namely: `): print(`1/I*(x-a1+b1*I)^(A*n1)*(x-a1-b1*I)^(A*n1)*(x-c1)^(A*n1)*(x-d1+b1*I)^(A*n1)*(x-d1-b1*I)^(A*n1)/(x^(B*n1+1)*(x-e1)^(B*n1+1)):`): print(`Try: `): print(`JF(4, 2, 5, 6, 10,3,5,n,x);`): elif nargs=1 and args[1]=JFSh then print(`JFSh(n1,x): the originral integrand of Salikohv. Try:`): print(`JSh(n1,x);`): elif nargs=1 and args[1]=JFShFe then print(`JShFe(n1):Fast version of JSh(2*n1), using the third-order linear recurrence with polynomial coefficients found via the`): print(`Almkvist-Zeilberger algorithm. Try:`): print(`JShFe(3);`): elif nargs=1 and args[1]=MakeDB then print(`MakeDB(M,n,N): inputs a positive integer M, and outputs a database of all the operators for computing`): print(`JF(4,2,5,6,10,A,B,2*n,x) for 1<=A,B<=M and gcd(A,B)=1. The output is a table such that T[A,B] is the desired operator.`): print(` Try: `): print(`MakeDB(2,n,N) ; `): elif nargs=1 and args[1]=JSh then print(`JSh(n1):The original integral, for odd n1 it gives combinations of 1 and arctan(1/7). For even n1 combinations of 1 and Pi. Try:`): print(`[seq(JSh(n1),n1=1..10)]`): elif nargs=1 and args[1]=NorPair then print(` NorPair(P,r,K): Given a pair of sequences of rational numbers P, multiplies the n-th term by `): print(` by lcm(1..r*n)*K^n. It also returns the list of gcds, let's call it G, and the list of log(G[i])/i and `): print(` and the smallest and largest log(G[i])/i in the second half, in floating point. Try: `): print(` NorPair(PairSeqPiAB(3,5,100),10,25/32); `): print(` NorPair(PairSeqPiAB(2,3,100),8,1/32); `): elif nargs=1 and args[1]=Nu1 then print(`Nu1(a,b,K,r): if you have a sequence of diophantine approximations A(n)+B(n)*c such that`): print(`(i)|A(n)+B(n)*c|<=Cb^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers`): print(`outputs the implied upper bound for the irrationality measure.`): print(`Here we do not use the fact that gcd(A(n),B(n)) may be of exponential growth.`): print(`For the more general case, try Nu1Extra(a,b,K,r,Extra); `): print(`Nu1(op(ShoreshPiAB(N,3,5)[2]),25/16,10);`): print(`Nu1(op(ShoreshPiAB(N,2,3)[2]), 1/32,8);`): elif nargs=1 and args[1]=Nu1Extra then print(`Nu1Extra(a,b,K,r,K1): if you have a sequence of diophantine approximations A(n)+B(n)*c such that`): print(`(i)|A(n)+B(n)*c|<=CONST*b^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers`): print(`AND it looks like, empirically or rigorously that gcd(K^n*dn(r*n)*A(n), K^n*dn(r*n)*B(n)) are O(e^(K1*n))`): print(`outputs the implied upper bound for the irrationality measure. Try:`): print(`Nu1Extra(op(ShoreshPiAB(N,3,5)[2]),25/16,10,0);`): print(`Nu1Extra(op(ShoreshPiAB(N,2,3)[2]),1/32,8,0);`): elif nargs=1 and args[1]=Ope23 then print(`Ope23(n,N): the linear recurrence equation annihilating the (2,3)-Pi-sequence, followed by the indicial equation in N,`): print(` and the roots. Try:`): print(`Ope23(n,N)`): elif nargs=1 and args[1]=OpePiAB then print(`OpePiAB(A,B,n,N): the linear difference operator in backward notation for computing JF(4,2,5,6,10,A,B,2*n,x). Try`): print(`OpePiAB(1,1,n,N); `): elif nargs=1 and args[1]=OpeSh then print(`OpeSh(n,N): the linear recurrence operator ope(n,N^(-1)) such that the Salikhov sequence c(n) is such that c(n)=ope(n,N^(-1))c(n).`): print(`it is for the full sequence, both odd and even. It is th-order. Try: `): print(`OpeSh(n,N);`): elif nargs=1 and args[1]=OpeShe then print(`OpeShe(n,N): the linear recurrence operator ope(n,N^(-1)) such that the Salikhov sequence c(2*n) is such that c(2*n)=ope(n,N^(-1))c(2*n).`): print(`Try: `): print(`OpeShe(n,N);`): elif nargs=1 and args[1]=OpeSho then print(`OpeSho(n,N): the linear recurrence operator ope(n,N^(-1)) such that the Salikhov sequence c(2*n) is such that c(2*n+1)=ope(n,N^(-1))c(2*n+1).`): print(`Try: `): print(`OpeShe(n,N);`): elif nops([args])=1 and op(1,[args])=PairSeqAT7 then print(`PairSeqAT7(N): the sequence of coefficients of 1 and arctan(1/7), respectively in AT7SeqF(N). Try:`): print(`PairSeqAT7(10);`): elif nargs=1 and args[1]=PairSeqPi then print(`PairSeqPi(N): the sequence of coefficients of 1 and Pi, respectively in PiSeqF(N). Try:`): print(`For a much faster version, use PairSeqPiAB(3,5,N); Try: `): print(`PairSeqPi(10);`): elif nargs=1 and args[1]=PairSeqPiAB then print(`PairSeqPiAB(A,B,N): the sequence of coefficients of 1 and Pi, respectively in PiSeqAB(N) (q.v.). Try:`): print(` VERY FAST, using pre-computed recurrences`): print(`for 1<=A,B<=10 and gcd(A,B)=1 and (A,B)=(8,13) ,. Otherwise it is the same speed.`): print(`PairSeqPiAB(3,5,100);`): elif nargs=1 and args[1]=Paper23 then print(`Paper23(M): inputs a large positive integer M, and outputs an article about diophantine approximations to Pi implied`): print(`by V. Kh. Salikhov integral (2) where (3,5) is replaced with (2,3).`): print(`"On the Measure of Irrationality of the Number Pi"`): print(`Mathematical Notes,2010, vol.88, no. 4, 563-573, but uses an alternative approach `): print(`using recurrences obtained by the Almkvist-Zeilberger algorithm. `): print(`It proves rigorously (modulo a simple divisibility lemma) a coarse upper bound for the irrationality measure `): print(`and a better one modulo a conjecture about the rate of growth of the gcd of the coefficients of 1 and Pi in the`): print(`relevant linear form. To get Salikhov's original theorem type. M is a large positive integer (at least 100) for the empirical determination.`): print(`Paper23(400); `): elif nargs=1 and args[1]=PaperAB then print(`PaperAB(A,B,M): inputs relatively prime positive integers A and B and a large positive integer M, at least 100`): print(`outputs an article about diophantine approximations to Pi implied`): print(`by V. Kh. Salikhov integral (2) in `): print(`"On the Measure of Irrationality of the Number Pi"`): print(`Mathematical Notes,2010, vol.88, no. 4, 563-573, but uses an alternative approach and where the exponets 3 and 5 in the`): print(`numerator and denominator are replaced by A and and B.`): print(`It proves rigorously (modulo a simple divisibility lemma) a coarse upper bound for the irrationality measure `): print(`and a better one modulo a conjecture about the rate of growth of the gcd of the coefficients of 1 and Pi in the`): print(`relevant linear form. To get Salikhov's original theorem type. M is a large positive integer (at least 100) for the empirical determination.`): print(`PaperAB(3,5,400); To get a conjecturally better bound, type:`): print(`PaperAB(2,3,400);`): elif nargs=1 and args[1]=PiSeq23 then print(`PiSeq23(M): The first M terms of the (2,3) analog of Salikhov's sequence, starting at n=1. Try:`): print(`PiSeq23(100):`): elif nargs=1 and args[1]=PiSeqF then print(`PiSeqF(N): the first N terms in the Salikhov even sequence, via the recurrence.`): print(`For an even faster, using the pre-computed recurrence (rather than finding the recurrence for scratch, use PiSeqAB(3,5,N);`): print(`Try: `): print(`PiSeqF(10);`): elif nargs=1 and args[1]=PiSeq then print(`PiSeq(M): the first M terms in the Salikhov even sequence. DONE DIRECTLY. VERY SLOW, only for checking. For much faster version use PiSeqF(M). Try:`): print(`PiSeq(3);`): elif nargs=1 and args[1]=PiSeqAB then print(`PiSeqAB(A,B,M): The first M terms of J(4,2,5,6,10,A,B,n0) for n0 from 1 to M. `): print(` VERY FAST, using pre-computed recurrences`): print(`for 1<=A,B<=10 and gcd(A,B)=1 and (A,B)=(8,13) ,. Otherwise it is the same speed.`): print(`Try:`): print(`PiSeqAB(3,5,30);`): elif nargs=1 and args[1]=PiSeqF then print(`PiSeqF(M): the first M terms in the Salikhov even sequence. DONE via the Almkvist-Zeilberger recurrence. Much faster. Try: `): print(`PiSeqF(20);`): elif nargs=1 and args[1]=SeqFromRec then print(`SeqFromRec(ope,n,N,INI,M): given a linear recurrence operator ope in n and the negative shift operator N^(-1) and`): print( ` (a list of initial values (starting at n=0), outputs the first M+1 terms of the sequence starting at n=0. `): pritn(` Try: `): print(` SeqFromRec(1/N+1/N^2,n,N,[0,1],10); `): elif nargs=1 and args[1]=ShoreshPiAB then print(`ShoreshPiAB(N,A,B): the indicial polynomial of the recurrence for the Salikhov Pi sequence`): print(`absolute values of the largest and smallest roots of the indicial equation of the recurrence for the Salikhov Pi sequence`): print(`linear recurrence operator with constant coefficients for Salikhov's approximation`): print(`Try:`): print(`ShoreshPiAB(N,2,3);`): print(`ShoreshPiAB(N,3,5);`): else print(`There is no ezra for`,args): fi: end: #########From EKHAD pashetd:=proc(p,N) local gu,p1,mekh,i: p1:=normal(p): mekh:=denom(p1): p1:=numer(p1): p1:=expand(p1): gu:=0: for i from 0 to degree(p1,N) do gu:=gu+factor(coeff(p1,N,i))*N^i: od: RETURN(gu,mekh): end: goremd:=proc(N,ORDER,bpc) local i, gu: gu:=0: for i from 0 to ORDER do gu:=gu+(bpc[i])*N^i: od: RETURN(gu): end: duisd:= proc(A,ORDER,y,n,N) local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu, mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1,bpc,apc: KAMA:=10: gorem:=goremd(N,ORDER,bpc): conj:=gorem: yakhas:=0: for i1 from 0 to degree(conj,N) do yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A): od: lu1a:=yakhas: LU1:=numer(yakhas): yakhas:=1/denom(yakhas): yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A): P:=LU1: Q:=numer(yakhas): R:=denom(yakhas): j1:=0: while j1 <= KAMA do g:=gcd(R,Q-j1*diff(R,y)): if g <> 1 then Q2:=(Q-j1*diff(R,y))/g: R:=normal(R/g): Q:=normal(Q2+j1*diff(R,y)): P:=P*g^j1: j1:=-1: fi: j1:=j1+1: od: P:=expand(P): R:=expand(R): Q:=expand(Q): Q1:=Q+diff(R,y): Q1:=expand(Q1): l1:=degree(R,y): l2:=degree(Q1,y): meka1:=coeff(R,y,l1): mekb1:=coeff(Q1,y,l2): if l1-1 <>l2 then k1:=degree(P,y)-max(l1-1,l2): else mumu:= -mekb1/meka1: if type(mumu,integer) and mumu > 0 then k1:=max(mumu, degree(P,y)-l2): else k1:=degree(P,y)-l2: fi: fi: fu:=0: if k1 < 0 then RETURN(0): fi: for ia1 from 0 to k1 do fu:=fu+apc[ia1]*y^ia1: od: gugu:=Q1*fu+R*diff(fu,y)-P: gugu:=expand(gugu): degg:=degree(gugu,y): for ia1 from 0 to degg do eq[ia1+1]:=coeff(gugu,y,ia1)=0: od: for ia1 from 0 to k1 do va1[ia1+1]:=apc[ia1]: od: for ia1 from 0 to ORDER do va2[ia1+1]:=bpc[ia1]: od: eq:=convert(eq,set): va1:=convert(va1,set): va2:=convert(va2,set): va1:=va1 union va2 : va1:=solve(eq,va1): fu:=subs(va1,fu): gorem:=subs(va1,gorem): if fu=0 and gorem=0 then RETURN(0): fi: for ia1 from 0 to k1 do gorem:=subs(apc[ia1]=1,gorem): od: for ia1 from 0 to ORDER do gorem:=subs(bpc[ia1]=1,gorem): od: fu:=normal(fu): shad:=pashetd(gorem,N): S:=lu1a*R*fu*shad[2]/P: S:=subs(va1,S): S:=normal(S): S:=factor(S): for ia1 from 0 to k1 do S:=subs(apc[ia1]=1,S): od: for ia1 from 0 to ORDER do S:=subs(bpc[ia1]=1,S): od: RETURN(shad[1],S): end: AZd:= proc(A,y,n,N) local ORDER,gu,KAMA: option remember: KAMA:=30: for ORDER from 0 to KAMA do gu:=duisd(A,ORDER,y,n,N): if gu<>0 then RETURN(gu): fi: od: 0: end: AZd1:= proc(A,y,n,N) local gu,ORDER,i1: option remember: gu:=AZd(A,y,n,N): if gu<>0 then gu:=gu[1]: ORDER:=degree(gu,N): gu:=gu/coeff(gu,N,ORDER): gu:=1-subs(n=n-ORDER,gu)/N^ORDER: gu:=add(factor(coeff(gu,N,-i1))*N^(-i1),i1=1..ORDER): RETURN(gu): fi: 0: end: #########end From EKHAD dn:=proc(n) local i: lcm(seq(i,i=1..n)): end: #deltE(a,c): Given a rational number a and a constant #c, finds the delta such that |a-c|=1/denom(a)^(1+delta) #For example, try delta(22/7,evalf(Pi)): deltE:=proc(a,c) local gu: if type(a,integer) then RETURN(FAIL): fi: gu:=evalf(-log(abs(c-a))/log(denom(a)))-1: evalf(gu,20): end: #JF(a1,b1,c1,d1,e1,A,B,n1,x) : the integrand of the generalized Salikhov integral. Try: #JF(4, 2, 5, 6, 10,3,5,n,x); JF:=proc(a1,b1,c1,d1,e1,A,B,n1,x) : 1/I*(x-a1+b1*I)^(A*n1)*(x-a1-b1*I)^(A*n1)*(x-c1)^(A*n1)*(x-d1+b1*I)^(A*n1)*(x-d1-b1*I)^(A*n1)/(x^(B*n1+1)*(x-e1)^(B*n1+1)): end: #J(a1,b1,c1,d1,e1,,A,B,n1): #1/I*int((x-a1+b1*I)^(A*n)*(x-a1-b1*I)^(A*n)*(x-c1)^(A*n)*(x-d1+b1*I)^(A*n)*(x-d1-b1*I)^(A*n)/(x^(B*n)*(x-e1)^(B*n)),x=a1-b1*I..a1+b1*I): #The original Salikhov integral is #J(4,2,5,6,10,3,5,n1); J:=proc(a1,b1,c1,d1,e1,A,B,n1) local x:simplify(int(JF(a1,b1,c1,d1,e1,A,B,n1,x),x=a1-b1*I..a1+b1*I)):end: #JFSh(n1,x): the originral integrand of Salikohv JFSh:=proc(n1,x): JF(4,2,5,6,10,3,5,n1,x): end: #JSh(n1):The original integral JSh:=proc(n1) option remember: J(4,2,5,6,10,3,5,n1): end: #OpeSh(n,N): the linear recurrence operator ope(n,N^(-1)) such that the Salikhov sequence c(n) is such that c(n)=ope(n,N^(-1))c(n). #Try: #OpeSh(n,N); OpeSh:=proc(n,N) local x: option remember: AZd1(JFSh(n,x),x,n,N): end: #OpeShe(n,N): the linear recurrence operator ope(n,N^(-1)) such that the Salikhov sequence c(2*n) is such that c(n)=ope(n,N^(-1))c(n). #Try: #OpeShe(n,N); OpeShe:=proc(n,N) local x: option remember: AZd1(JFSh(2*n,x),x,n,N): end: #ShoreshPi(N): the indicial polynomial of the recurrence for the Salikhov Pi sequence #absolute values of the largest and smallest roots of the indicial equation of the recurrence for the Salikhov Pi sequence #linear recurrence operator with constant coefficients for Salikhov's approximation #Try: #ShoreshPi(N); ShoreshPi:=proc(N) local x,ope,lu,i,n: option remember: ope:=AZd(JFSh(2*n,x),x,n,N)[1]: ope:=lcoeff(ope,n): ope:=normal(ope/igcd(seq(coeff(ope,N,i),i=0..degree(ope,N)))): lu:=[solve(ope,N)]: lu:=[seq(abs(lu[i]),i=1..nops(lu))]: ope, evalf([max(abs(lu[1]),abs(lu[2]) ), min(abs(lu[1]),abs(lu[2]))],20): end: #OpeSho(n,N): the linear recurrence operator ope(n,N^(-1)) such that the Salikhov sequence c(2*n+1) is such that c(n)=ope(n,N^(-1))c(n). #Try: #OpeSho(n,N); OpeSho:=proc(n,N) local x: option remember: AZd1(JFSh(2*n+1,x),x,n,N): end: ###########Start Pi #PiSeq(N): the first N terms in the Salikhov even sequence. Try: #PiSeq(3); PiSeq:=proc(N) local n1: [seq(JSh(2*n1),n1=1..N)]: end: #JShFe(n1):Fast version of JSh(2*n1), using the third-order linear recurrence with polynomial coefficients found via the #Almkvist-Zeilberger algorithm. Try: #JShFe(3); JShFe:=proc(n1) local ope,i1,n,N: option remember: if not (type(n1,integer) and n1>=0) then print(n1, `must be a positive integer`): RETURN(FAIL): fi: if n1<=2 then RETURN(JSh(2*n1)): fi: ope:=OpeShe(n,N): add(subs(n=n1,coeff(ope,N,-i1))*JShFe(n1-i1),i1=1..-ldegree(ope,N)): end: #PiSeqF(N): the first N terms in the Salikhov even sequence. Try: #PiSeqF(3); PiSeqF:=proc(N) local n1: option remember: [seq(JShFe(n1),n1=1..N)]: end: #PairSeqPi(N): the sequence of coefficients of 1 and Pi, respectively in PiSeqF(N). Try: #PairSeqPi(10); PairSeqPi:=proc(N) local gu,i: option remember: gu:=PiSeqF(N): [[seq(coeff(gu[i],Pi,0),i=1..N)],[seq(coeff(gu[i],Pi,1),i=1..N)]]: end: #AppSeqPi(N): the first N terms of the Salikhov approximating sequence to Pi. Try: #AppSeqPi(20); AppSeqPi:=proc(N) local gu,i: option remember: gu:=PairSeqPi(N): [seq(-gu[1][i]/gu[2][i],i=1..N)]: end: #AppSeqPiAB(A,B,N): the first N terms of the generalized Salikhov approximating sequence to Pi using SeqPiAB(A,B,N). Try: #AppSeqPABi(3,5,100); AppSeqPiAB:=proc(A,B,N) local gu,i: option remember: gu:=PairSeqPiAB(A,B,N): [seq(-gu[1][i]/gu[2][i],i=1..N)]: end: ###########end Pi ###########Start arctan(1/7) #AT7Seq(N): the first N terms in the Salikhov odd sequence, featuring arctan(1/7). Try: #AT7Seq(10); AT7Seq:=proc(N) local n1: [seq(JSh(2*n1+1),n1=0..N)]: end: #JShFo(n1):Fast version of JSh(2*n1+1), using the third-order linear recurrence with polynomial coefficients found via the #Almkvist-Zeilberger algorithm. Try: #JShFo(3); JShFo:=proc(n1) local ope,i1,n,N: option remember: if not (type(n1,integer) and n1>=0) then print(n1, `must be a positive integer`): RETURN(FAIL): fi: if n1<=2 then RETURN(JSh(2*n1+1)): fi: ope:=OpeSho(n,N): add(subs(n=n1,coeff(ope,N,-i1))*JShFo(n1-i1),i1=1..-ldegree(ope,N)): end: #AT7SeqF(N): the first N terms in the Salikhov odd sequence, featuring arctan(1/7). Try: #AT7SeqF(3); AT7SeqF:=proc(N) local n1: option remember: [seq(JShFo(n1),n1=0..N)]: end: #PairSeqAT7(N): the sequence of coefficients of 1 and arctan(1/7), respectively in AT7SeqF(N). Try: #PairSeqAT7(10); PairSeqAT7:=proc(N) local gu,i: option remember: gu:=AT7SeqF(N): [[seq(coeff(gu[i], arctan(1/7) ,0),i=1..N)],[seq(coeff(gu[i],arctan(1/7),1),i=1..N)]]: end: #AppSeqAT7(N): the first N terms of the approximating sequence to arctan(1/7), using the Salikhov sequence. Try: #AppSeqAT7(20); AppSeqAT7:=proc(N) local gu,i: option remember: gu:=PairSeqAT7(N): [seq(-gu[1][i]/gu[2][i],i=1..N)]: end: ###########end arctan(1/7) #JFu(a1,b1,c1,d1,e1,A,B,n1,u) : the integrand of the generalized Salikhov integral when transferred to (0,1) #JFu(4, 2, 5, 6, 10,3,5,n,u); JFu:=proc(a1,b1,c1,d1,e1,A,B,n1,u) local gu,x: gu:=1/I*(x-a1+b1*I)^(A*n1)*(x-a1-b1*I)^(A*n1)*(x-c1)^(A*n1)*(x-d1+b1*I)^(A*n1)*(x-d1-b1*I)^(A*n1)/(x^(B*n1+1)*(x-e1)^(B*n1+1)): gu:=simplify(2*b1*I*subs(x=a1-b1*I+2*b1*I*u,gu)): end: #GuessK1(L,m): Inputs a sequence of integers L, # finds a constant K of the form K=ithprime(i)^j, with i and j<=m such that L[n]/K^n is an integer for all n. Try: #GuessK1([seq((i^3+1)*25^(i),i=1..30)],10); GuessK1:=proc(L,m) local K,i,lu,i1,p,K1,j: K:=1: for i from 1 to m do p:=ithprime(i): for j from m by -1 while max(seq(denom(L[i1]/(p^j)^i1),i1=1..nops(L)))>3000 to 0 do od: K:=K*p^j : od: K: end: #GuessK(L,m,r): inputs a sequence of rational numbers, L, a positive integer m, and a positive integer r #finds a rational number K such that L[i]*K^i*dn(r*i) are all integers. Try: #GuessK([seq((2/3)^i/dn(5*i),i=1..20)],30,5); GuessK:=proc(L,m,r) local L1,L1t,L1b,i,K1,K2,K,vu: L1:=[]: for i from 1 to nops(L) do L1:=[op(L1), L[i]*dn(r*i)]: od: L1t:=numer(L1): L1b:=denom(L1): K1:=GuessK1(L1t,m): K2:=GuessK1(L1b,m): K:=K2/K1: vu:= max(seq(denom(L[i]*dn(r*i)*K^i),i=1..nops(L)) ): if vu>3000 then #print(K, `did not work out`): RETURN(FAIL): fi: K: end: #Nu1(a,b,K,r): if you have a sequence of diophantine approximations A(n)+B(n)*c such that #(i)|A(n)+B(n)*c|<=Cb^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers #outputs the implied upper bound for the irrationality measure. Try: #Nu1(op(ShoreshPi(N)[2]),25/16,10); Nu1:=proc(a,b,K,r) local lu: lu:=-(log(b)+r+log(K))/(log(a)+r+log(K)): evalf(1+1/lu,20): end: #OpePiAB(A,B,n,N): the linear difference operator in backward notation for computing JF(4,2,5,6,10,A,B,2*n,x) OpePiAB:=proc(A,B,n,N) local x: if A<=10 and B<=10 then DB(A,B,n,N): elif A=8 and B=13 then Ope813(n,N): else AZd1(JF(4,2,5,6,10,A,B,2*n,x),x,n,N): fi: end: #MakeDB(M,n,N): inputs a positive integer M, and outputs a database of all the operators for computing #JF(4,2,5,6,10,A,B,2*n,x) for 1<=A,B<=N and gcd(A,B)=1. The output is a table such that T[A,B] is the desired operator. #Try: #MakeDB(2,n,N) MakeDB:=proc(M,n,N) local T,A,B: for A from 1 to M do for B from 1 to M do if igcd(A,B)=1 then T[A,B]:=OpePiAB(A,B,n,N): fi: od: od: op(T): end: #SeqFromRec(ope,n,N,INI,M): given a linear recurrence operator ope in n and the negative shift operator N^(-1) and #a list of initial values (starting at n=0), outputs the first M+1 terms of the sequence starting at n=0. #Try: #SeqFromRec(1/N+1/N^2,n,N,[0,1],10); SeqFromRec:=proc(ope,n,N,INI,M) local T,n0,d,i: d:=-ldegree(ope,N): if nops(INI)<>d then RETURN(FAIL): fi: for n0 from 0 to d-1 do T[n0]:=INI[n0+1]: od: for n0 from d to M do T[n0]:=add(subs(n=n0,coeff(ope,N,-i))*T[n0-i],i=1..d): od: [seq(T[n0],n0=0..M)]: end: #PiSeqAB(A,B,M): The first M terms of J(4,2,5,6,10,A,B,n0) for n0 from 1 to M. VERY FAST, using pre-computed recurrences #for 1<=A,B<=10 and gcd(A,B)=1 and (A,B)=(8,13) ,. Otherwise it is the same speed. #Try: #PiSeqAB(1,1,30); PiSeqAB:=proc(A,B,M) local INI, n0,ope,n,N: ope:=OpePiAB(A,B,n,N): INI:=[seq( J(4,2,5,6,10,A,B,2*n0),n0=0..2)]: [op(2..M+1,SeqFromRec(ope,n,N,INI,M))]: end: #EmpDelAB(A,B,M): the empirical deltas from M-10 to M for the Diophatine approximations implied by #J(4,2,5,6,10,A,B,2*n0). Try: #EmpDelAB(1,1,40); EmpDelAB:=proc(A,B,M) local gu,lu,i: gu:=PiSeqAB(A,B,M): lu:=[seq(deltE(-coeff(gu[i],Pi,0)/coeff(gu[i],Pi,1),Pi), i=M-10..M)]: evalf(lu,20): end: #ShoreshPiAB(N,A,B): the indicial polynomial of the recurrence for the Salikhov Pi sequence #absolute values of the largest and smalleat roots of the indicial equation of the recurrence for the Salikhov Pi sequence #linear recurrence operator with constant coefficients for Salikhov's approximation #Try: #ShoreshPiAB(N,2,3); #ShoreshPiAB(N,3,5); ShoreshPiAB:=proc(N,A,B) local x,ope,lu,i,n: option remember: ope:=AZd(JF(4,2,5,6,10,A,B,2*n,x),x,n,N)[1]: ope:=lcoeff(ope,n): ope:=normal(ope/igcd(seq(coeff(ope,N,i),i=0..degree(ope,N)))): lu:=[solve(ope,N)]: lu:=[seq(abs(lu[i]),i=1..nops(lu))]: ope, evalf([max(abs(lu[1]),abs(lu[2]) ), min(abs(lu[1]),abs(lu[2]))],20): end: #Nu1(a,b,K,r): if you have a sequence of diophantine approximations A(n)+B(n)*c such that #(i)|A(n)+B(n)*c|<=Cb^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers #outputs the implied upper bound for the irrationality measure. Try: #Nu1(op(ShoreshPi(N)[2]),25/16,10); Nu1:=proc(a,b,K,r) local lu: lu:=-(log(b)+r+log(K))/(log(a)+r+log(K)): evalf(1+1/lu,20): end: #PairSeqPiAB(A,B,N): the sequence of coefficients of 1 and Pi, respectively in PiSeqAB(N) (q.v.). Try: #PairSeqPiAB(3,5,10); PairSeqPiAB:=proc(A,B,N) local gu,i: option remember: gu:=PiSeqAB(A,B,N): [[seq(coeff(gu[i],Pi,0),i=1..N)],[seq(coeff(gu[i],Pi,1),i=1..N)]]: end: #GuessKr(L,m): inputs a sequence of rational numbers, L, a positive integer m, #finds a positive integer r, and a rational number K such that L[i]*K^i*dn(r*i) are all integers. Try: #GuessKr([seq((2/3)^i/dn(5*i),i=1..20)],30); GuessKr:=proc(L,m) local r,i,K: if nops(L)<20 then print(`List too short, make it at least of length 20`): RETURN(FAIL): fi: if type(L[nops(L)],integer) then r:=0: else r:=round(ifactors(denom(L[nops(L)]))[2][-1][1]/nops(L)); for i from 15 to nops(L) do if denom(L[i])<>1 then if not member(round(ifactors(denom(L[i]))[2][-1][1]/i),{r,r-1}) then # print(r, `did not work out`): RETURN(FAIL): fi: fi: od: fi: K:=GuessK(L,m,r): if K=FAIL then # print(`r is`, r, `but could not find K`): RETURN(FAIL): fi: [r,K]: end: #GuessKrPair(P,m): Given a pair,P, of sequences of rational numbers and a positive integer m #finds a non-negative integer r and a rational number K such that #BOTH P[1][n]*dn(r*n)*K^(i*n) and P[1][n]*dn(r*n)*K^(i*n) are sequences of INTEGERS. #It tries out all the primes up to m and all the powers up to m. #Type: GuessKrPair:=proc(P,m) local gu1,gu2: gu1:=GuessKr(P[1],m): if gu1=FAIL then RETURN(FAIL): fi: gu2:=GuessKr(P[2],m): if gu2=FAIL then RETURN(FAIL): fi: [max(gu1[1],gu2[1]), lcm(numer(gu1[2]), numer(gu2[2]))/igcd(denom(gu1[2]), denom(gu2[2]))]: end: #NorPair(P,r,K): Given a pair of sequences of rational numbers P, multiplies the n-th term by #by lcm(1..r*n)*K^n. It also returns the list of gcds, let's call it G, and the list of log(G[i])/i and #and the smallest and largest log(G[i])/i in the second half, in floating point. Try: #NorPair(PairSeqPiAB(3,5,100),10,25/32); NorPair:=proc(P,r,K) local P1,P2,vu,gu,i1: if not (type(P,list) and nops(P)=2 and type(P[1],list) and type(P[2],list) and nops(P[1])=nops(P[2])) then print(`Bad input`): RETURN(FAIL): fi: if not (type(r,integer) and r>=0 and (type(K,integer) or type(K,fraction)) ) then print(`Bad input`): RETURN(FAIL): fi: P1:=[seq(P[1][i1]*dn(r*i1)*K^i1,i1=1..nops(P[1]))]: P2:=[seq(P[2][i1]*dn(r*i1)*K^i1,i1=1..nops(P[2]))]: vu:=max(op(denom(P1)),op(denom(P2))): if vu>3000 then RETURN(FAIL): fi: P1:=vu*P1: P2:=vu*P2: gu:=[seq(igcd(P1[i1],P2[i1]),i1=1..nops(P1))]: vu:=[seq(evalf(log(gu[i1])/i1, 20),i1=trunc(nops(P1)/2)..nops(P1))]: [[P1,P2],gu,vu,min(op(vu)),max(op(vu))]: end: #Nu1Extra(a,b,K,r,K1): if you have a sequence of diophantine approximations A(n)+B(n)*c such that #(i)|A(n)+B(n)*c|<=CONST*b^n, (ii) max(|A(n)|,|B(n)|)=OMEGA(a^n) (iii) K^n*dn(r*n)*A(n) and K^n*dn(r*n)*B(n) are integers #AND it looks like, empirically or rigorously that gcd(K^n*dn(r*n)*A(n), K^n*dn(r*n)*B(n)) are O(e^(K1*n)) #outputs the implied upper bound for the irrationality measure. Try: #Nu1Extra(op(ShoreshPi(N)[2]),25/16,10,0); Nu1Extra:=proc(a,b,K,r,K1) local lu: lu:=-(log(b)+r+log(K)-K1)/(log(a)+r+log(K)-K1): evalf(1+1/lu,20): end: #InfoAB(A,B,x,n,N,M): inputs positive integers A and B such that gcd(A,B), symbols x, n, N and a large positive integer M #and outputs a list consisting of #(i)1/I*Int((x-4+2*I)^(2*A*n)*(x-4-2*I)^(2*A*n)*(x-5)^(2*A*n)*(x-6+2*I)^(2*A*n)*(x-6-2*I)^(2*A*n)/(x^(2*B*n+1)*(x-10)^(2*B*n+1)),x=4-2*I..4+2*I): #(this is the generalized Salikhov integral) that can be written as C(n)+D(n)*Pi for some rational numbers C(n) and D(n)/ #(ii) the Empirical delta for the approximation to Pi from the sequence obtained from the generalized Salikhov integral #with A, B i.e. #(iii) the linear recurrence operator in n and the shift operator N annihilating E(n) and hence the coefficients of Pi and 1 #(iv) The leading coefficient in n, in other words the constant-coefficient operator approximinating it. #the so-called indicial equation #(v) The pair [largest root, absolute value of smaller roots] #(vi) The conjectured pair [r,K] such that lcm(1..r*n)*K^n*E(n) has integer coefficients #Let's call the coefficients if 1 and Pi A1(n) and A2(n) respectively/ #(vii): The upper bound for the irrationality NOT taking advantage of the fact that gcd(C(n),D(n)) may #(vii) the smallest log(gcd(A1(n),A2(n))/n for n between M/2 and M #(viii) assuming that this is a lower bound for the lim sup of log(gcd(A1(n),A2(n))/n for n between M/2 and M #the implied upper bound for the irrationality measure of Pi. Try: #InfoAB(3,5,n,N,400); InfoAB:=proc(A,B,x,n,N,M) local ope,ope1,gu,lu0,lu1,lu2,lu3,lu4,lu5: option remember: if not (type(A,integer) and type(B,integer) and igcd(A,B)=1 and type(n,symbol) and type(N,symbol) and type(M,integer)) then print(`Bad input`): RETURN(FAIL): fi: gu:=[1/I*Int((x-4+2*I)^(2*A*n)*(x-4-2*I)^(2*A*n)*(x-5)^(2*A*n)*(x-6+2*I)^(2*A*n)*(x-6-2*I)^(2*A*n)/(x^(2*B*n+1)*(x-10)^(2*B*n+1)),x=4-2*I..4+2*I)]: lu0:=min(op(EmpDelAB(A,B,M))): lu0:=evalf(1+1/lu0,20): gu:=[op(gu),lu0]: ope:=AZd(JF(4,2,5,6,10,A,B,2*n,x),x,n,N)[1]: ope1:=lcoeff(ope,n): ope1:=numer(ope1/coeff(ope1,N,0)): gu:=[op(gu),ope,ope1]: lu1:=ShoreshPiAB(N,A,B)[2]: gu:=[op(gu),lu1]: lu2:=GuessKrPair(PairSeqPiAB(A,B,40),20); if lu2=FAIL then RETURN(gu): else gu:=[op(gu),lu2]: fi: lu3:=Nu1Extra(lu1[1],lu1[2],lu2[2],lu2[1],0): gu:=[op(gu),lu3]: lu4:=NorPair(PairSeqPiAB(A,B,M),lu2[1],lu2[2])[4]; gu:=[op(gu),lu4]: lu5:=Nu1Extra(lu1[1],lu1[2],lu2[2],lu2[1],lu4): gu:=[op(gu),lu5]: gu: end: #PaperAB(A,B,M): inputs relatively prime positive integers A and B and outputs an article about diophantine approximations to Pi implied #by V. Kh. Salikhov integral (2) in #"On the Measure of Irrationality of the Number Pi"`): #Mathematical Notes,2010, vol.88, no. 4, 563-573, but uses an alternative approach and where the exponets 3 and 5 in the #numerator and denominator are replaced by A and and B. #It proves rigorously (modulo a simple divisibility lemma) a coarse upper bound for the irrationality measure #and a better one modulo a conjecture about the rate of growth of the gcd of the coefficients of 1 and Pi in the #relevant linear form. To get Salikhov's original theorem type. M is a large positive integer (at least 100) for the empirical determination. #PaperAB(3,5,400); To get a conjecturally better bound, type: #PaperAB(2,3,400); PaperAB:=proc(A,B,M) local x,n,N,gu,ope,ope1,F,i1,a,b,C,D,C1,D1,r,K,a1,b1,lu,K1,t0,n0: t0:=time(): gu:=InfoAB(A,B,x,n,N,M): if gu=FAIL or nops(gu)<9 then print(`Sorry, no paper`): RETURN(FAIL): fi: if gu[9]>0 then print(`A sequence of Diophantine Approximations to Pi that implies its irrationality`): print(`With a fully rigorous crude upper bound on the irrationality measure of`, gu[7] ): print(`and a not-yet-rigorous smaller upper bound around`, gu[9]): print(``): print(`By Shalosh B. Ekhad `): print(``): else print(`A sequence of Diophantine Approximations to Pi that Disappointingly does not imply its irrationality`): print(``): print(`By Shalosh B. Ekhad `): print(``): fi: print(`Consider the following Salikhov-type integral, let's call it E(n) `): print(``): print(gu[1]): print(``): print(`It is readily seen that E(n) can be written as`): print(``): print(`E(n)=C(n) + D(n)*Pi `): print(``): print(`For certain sequences of RATIONAL numbers C(n), D(n)`): print(``): ope:=gu[3]: print(`It follows from the amazing Almkvist-Zeilberger algorithm that E(n), and hence C(n) and D(n) satisfy the following`): print(` third-order linear recurrence equation with polynomial coefficients. `): print(``): print(add(coeff(ope,N,i1)*F(n+i1),i1=0..degree(ope,N))=0): print(``): print(`and in Maple format `): print(``): lprint(add(coeff(ope,N,i1)*F(n+i1),i1=0..degree(ope,N))=0): print(``): print(`The initial conditions are as follows`): print(``): print(seq( E(n0)=J(4,2,5,6,10,A,B,2*n0),n0=0..2)): print(``): print(`and in Maple notation`): print(``): lprint(seq( E(n0)=J(4,2,5,6,10,A,B,2*n0),n0=0..2)): print(``): print(`Using this recurrence, we can compute many values, and use them to ESTIMATE the implied bound irrationality measure.`): print(`Using the values from n=`, M-10, ` to M=`, M, `and taking the largest value yields the following estimate`): print(``): lprint(gu[2]): print(``): print(`We can use the Poincare lemma to find the asymptotic behavior of the exponentially growing C(n), D(n) and`): print(`the exponentially decaying E(n).`): print(``): ope1:=gu[4]: print(`The constant-coefficient approximation to the above recurrence has indicial polynomial `): print(``): print(ope1): print(``): print(`and in Maple format `): print(``): lprint(ope1): print(``): print(`whose largest root let's call it a, is`, gu[5][1], `and absolute value of the two smaller roots, let's call it b is`, gu[5][2]): print(``): print(`It follows that, up to polynomial correction that will get out in the wash, C(n) and D(n) are OMEGA(a^n) and E(n) is O(b^n)`): print(``): print(`We now need a divisibility lemma that is left to the reader`): print(``): print(`Let d(n)=lcm(1..n) `): r:=gu[6][1]: K:=gu[6][2]: print(`Lemma: `): print(C1(n)=d(r*n)*K^n*C(n), D1(n)=d(r*n)*K^n*D(n)): print(`are integers`): print(`Defining `): print(``): print(E1(n)=d(r*n)*K^n*E(n)): print(``): print(`We have `): print(``): print(E1(n)=C1(n)+ Pi*D1(n)): print(``): print(`where C1(n) and D1(n) are INTEGERS `): print(``): print(`Since, famously, d(n)=OMEGA(e^n), we have `): print(``): print(`C1(n) and D1(n) are OMEGA of`): print(exp(r*n)*K^n*a^n): print(``): print(`and E1(n) is O of `): print(``): print(exp(r*n)*K^n*b^n): print(``): a1:=gu[5][1]: b1:=gu[5][2]: print(`It follows that E1(n)=max(C1(n),D1(n))^(-delta) where delta equals `): print(``): lu:=-(log(b)+r+log(K))/(log(a)+r+log(K)): print(lu): print(`recall that a is `, a1, `and b is `, b1): lu:=evalf(subs({b=b1,a=a1},lu),20): print(`and delta is `): print(``): lprint(lu): print(``): print(`implying the rigorous, but crude upper bound for the irrationality measure, 1+1/delta that equals `): print(``): lprint(evalf(1+1/lu,20)): print(``): print(`confirming the first part of the statement in the title `): print(``): print(`We now need a harder lemma, that we confess that we can't do, but you the reader may be able to `): print(``): print(`Harder Lemma: There exists a rigorously proved constant K1 such that `): print(``): print(`gcd(C1(n),D1(n))=O(exp(K1*n)) `): print(``): print(`Once we find such a constant K1 the delta can be improved to `): print(``): lu:=-(log(b)+r+log(K)-K1)/(log(a)+r+log(K)-K1): print(``): print(`By looking at the smallest value of log(gcd(C1(n),D1(n))/n for n between n=`, M/2, `and `, M, `we estimate that K1 can be taken to be`): print(``): lprint(gu[8]): print(``): print(`and the better delta is `): print(``): lu:=subs({a=a1,b=b1,K1=gu[8]},lu): lu:=evalf(lu,20): print(``): lprint(lu): print(``): print(`that implies an irrationality measure, 1+1/delta that equals `): print(``): lprint(evalf(1+1/lu,20)): print(``): print(`-----------------------------------------`): print(``): print(`This ends this paper that took`, time()-t0, `to generate. `): print(``): end: #Ope23(n,N): the linear recurrence equation annihilating the (2,3)-Pi-sequence, followed by the indicial equation in N, #and the roots. Try: #Ope23(n,N) Ope23:=proc(M) local ope,ope1,gu: ope:=AZd(1/I*(x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x,n,N)[1]: ope1:=lcoeff(ope,n): ope1:=numer(ope1/coeff(ope1,N,0)): gu:=evalf([solve(ope1,N)],20): [ope,ope1,gu]: end: #PiSeq23(M): The first M terms of the (2,3) analog of Salikhov's sequence, starting at n=1. Try: #PiSeq23(100): PiSeq23:=proc(M) local INI,ope,n,N,n1,gu: INI:=[ seq(1/I*int((x-4+2*I)^(2*n1)*(x-4-2*I)^(2*n1)*(x-5)^(2*n1)*(x-6+2*I)^(2*n1)*(x-6-2*I)^(2*n1)/(x^(3*n1+1)*(x-10)^(3*n1+1)),x=4-2*I..4+2*I), n1=0..2)]: INI:=simplify(INI): ope:=AZd1(1/I*(x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x,n,N): gu:=SeqFromRec(ope,n,N,INI,M): [op(2..nops(gu),gu)]: end: #Info23(x,n,N,M): inputs symbols x,n,N, and a large positive integer N, and outputs a list consisting of #(i)1/I*Int((x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*A*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*A*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x=4-2*I..4+2*I): #(this is the generalized Salikhov integral) that can be written as C(n)+D(n)*Pi for some rational numbers C(n) and D(n)/ #(ii) the Empirical delta for the approximation to Pi from the sequence obtained from the modified Salikhov integral #with 3 replaced by 2 and 5 replaced by 3, but vaild for all n, not just even n. #(iii) the linear recurrence operator in n and the shift operator N annihilating E(n) and hence the coefficients of Pi and 1 #(iv) The leading coefficient in n, in other words the constant-coefficient operator approximinating it. #the so-called indicial equation #(v) The pair [largest root, absolute value of smaller roots] #(vi) The conjectured pair [r,K] such that lcm(1..r*n)*K^n*E(n) has integer coefficients #Let's call the coefficients if 1 and Pi A1(n) and A2(n) respectively/ #(vii): The upper bound for the irrationality NOT taking advantage of the fact that gcd(C(n),D(n)) may #(vii) the smallest log(gcd(A1(n),A2(n))/n for n between M/2 and M #(viii) assuming that this is a lower bound for the lim sup of log(gcd(A1(n),A2(n))/n for n between M/2 and M #the implied upper bound for the irrationality measure of Pi. Try: #Info23(x,n,N,400); Info23:=proc(x,n,N,M) local ope,opeNeg,INI,Sid, Sid1, ope1,gu,lu0,lu1,lu2,lu3,lu4,lu5,SidA,SidB,vu: option remember: if not ( type(n,symbol) and type(N,symbol) and type(M,integer)) then print(`Bad input`): RETURN(FAIL): fi: gu:=[1/I*Int((x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x=4-2*I..4+2*I)]: ope:=AZd(1/I*(x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x,n,N)[1]: opeNeg:=AZd1(1/I*(x-4+2*I)^(2*n)*(x-4-2*I)^(2*n)*(x-5)^(2*n)*(x-6+2*I)^(2*n)*(x-6-2*I)^(2*n)/(x^(3*n+1)*(x-10)^(3*n+1)),x,n,N): INI:=[ seq(1/I*int((x-4+2*I)^(2*n1)*(x-4-2*I)^(2*n1)*(x-5)^(2*n1)*(x-6+2*I)^(2*n1)*(x-6-2*I)^(2*n1)/(x^(3*n1+1)*(x-10)^(3*n1+1)),x=4-2*I..4+2*I), n1=0..2)]: INI:=simplify(INI): Sid:=SeqFromRec(opeNeg,n,N,INI,M): Sid1:=[op(nops(Sid)-20..nops(Sid),Sid)]: Sid1:=[seq(deltE(-coeff(Sid1[i],Pi,0)/coeff(Sid1[i],Pi,1),Pi),i=1..nops(Sid1))]: lu0:=min(op(Sid1)): lu0:=evalf(1+1/lu0,20): gu:=[op(gu),lu0,ope]: ope1:=lcoeff(ope,n): ope1:=numer(ope1/coeff(ope1,N,0)): gu:=[op(gu),ope1]: lu1:=[solve(ope1,N)]: lu1:=evalf(lu1): lu1:=[lu1[1],abs(lu1[2])]: lu1:=evalf(lu1,20): gu:=[op(gu),lu1]: vu:=lcm(seq(denom(Sid[i]*dn(4*i)/32^trunc(i/2)),i=1..nops(Sid))): if vu>100 then RETURN(FAIL): else gu:=[op(gu),[4,1/(4*sqrt(2))]]: fi: lu3:=Nu1Extra(lu1[1],lu1[2],1/(4*sqrt(2)),4,0): gu:=[op(gu),lu3]: SidA:=[seq(coeff(vu*Sid[i],Pi,0),i=1..nops(Sid))]: SidB:=[seq(coeff(vu*Sid[i],Pi,1),i=1..nops(Sid))]: SidA:=[seq(dn(4*i)/32^trunc(i/2)*SidA[i],i=1..nops(SidA))]: SidB:=[seq(dn(4*i)/32^trunc(i/2)*SidB[i],i=1..nops(SidB))]: if max(op(denom(SidA)),op(denom(SidB)))>1000 then print(`Something terrible happened`): RETURN(FAIL): fi: lu4:=[seq(log(igcd(SidA[i],SidB[i]))/i,i=trunc(M/2)..M)]: lu4:=min(lu4): lu4:=evalf(lu4,20): gu:=[op(gu),lu4]: lu5:=Nu1Extra(lu1[1],lu1[2],1/(4*sqrt(2)),4,lu4): gu:=[op(gu),lu5]: gu: end: #Paper23(M): inputs a large positive integer M, and outputs an article about diophantine approximations to Pi implied #by V. Kh. Salikhov integral (2) where (3,5) is replaced with (2,3). #"On the Measure of Irrationality of the Number Pi"`): #Mathematical Notes,2010, vol.88, no. 4, 563-573, but uses an alternative approach #It proves rigorously (modulo a simple divisibility lemma) a coarse upper bound for the irrationality measure #and a better one modulo a conjecture about the rate of growth of the gcd of the coefficients of 1 and Pi in the #relevant linear form. To get Salikhov's original theorem type. M is a large positive integer (at least 100) for the empirical determination. #Paper23(400); Paper23:=proc(M) local x,n,N,gu,ope,ope1,F,i1,a,b,C,D,C1,D1,a1,b1,lu,K1,t0,n0: t0:=time(): gu:=Info23(x,n,N,M): print(`A sequence of Diophantine Approximations to Pi that implies its irrationality`): print(`With a fully rigorous crude upper bound on the irrationality measure of`, gu[7] ): print(`and a not-yet-rigorous smaller upper bound around`, gu[9]): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Consider the following Salikhov-type integral, let's call it E(n) `): print(``): print(gu[1]): print(``): print(`It is readily seen that E(n) can be written as`): print(``): print(`E(n)=C(n) + D(n)*Pi `): print(``): print(`For certain sequences of RATIONAL numbers C(n), D(n)`): print(``): ope:=gu[3]: print(`It follows from the amazing Almkvist-Zeilberger algorithm that E(n), and hence C(n) and D(n) satisfy the following`): print(` third-order linear recurrence equation with polynomial coefficients. `): print(``): print(add(coeff(ope,N,i1)*F(n+i1),i1=0..degree(ope,N))=0): print(``): print(`and in Maple format `): print(``): lprint(add(coeff(ope,N,i1)*F(n+i1),i1=0..degree(ope,N))=0): print(``): print(`The initial conditions are as follows`): print(``): print(seq( E(n0)=J(4,2,5,6,10,2,3,n0),n0=0..2)): print(``): print(`and in Maple notation`): print(``): lprint(seq( E(n0)=J(4,2,5,6,10,2,3,n0),n0=0..2)): print(``): print(`Using this recurrence, we can compute many values, and use them to ESTIMATE the implied bound irrationality measure.`): print(`Using the values from n=`, M-10, ` to M=`, M, `and taking the largest value yields the following estimate`): print(``): lprint(gu[2]): print(``): print(`We can use the Poincare lemma to find the asymptotic behavior of the exponentially growing C(n), D(n) and`): print(`the exponentially decaying E(n).`): print(``): ope1:=gu[4]: print(`The constant-coefficient approximation to the above recurrence has indicial polynomial `): print(``): print(ope1): print(``): print(`and in Maple format `): print(``): lprint(ope1): print(``): print(`whose largest root let's call it a, is`, gu[5][1], `and absolute value of the two smaller roots, let's call it b is`, gu[5][2]): print(``): print(`It follows that, up to polynomial correction that will get out in the wash, C(n) and D(n) are OMEGA(a^n) and E(n) is O(b^n)`): print(``): print(`We now need a divisibility lemma that is left to the reader`): print(``): print(`Let d(n)=lcm(1..n) `): print(`Lemma: `): print(C1(n)=d(4*n)/32^(trunc(n/2))*C(n), D1(n)=d(4*n)/32^(trunc(n/2))*D(n)): print(`are integers`): print(`Defining `): print(``): print(E1(n)=d(4*n)/(4*sqrt(2))^n*E(n)): print(``): print(`We have `): print(``): print(E1(n)=C1(n)+ Pi*D1(n)): print(``): print(`where C1(n) and D1(n) are INTEGERS `): print(``): print(`Since, famously, d(n)=OMEGA(e^n), we have `): print(``): print(`C1(n) and D1(n) are OMEGA of`): print(exp(4*n)/(4*sqrt(2))^n*a^n): print(``): print(`and E1(n) is Big O of `): print(``): print(exp(4*n)/(4*sqrt(2))^n*b^n): print(``): a1:=gu[5][1]: b1:=gu[5][2]: print(`It follows that E1(n)=max(C1(n),D1(n))^(-delta) where delta equals `): print(``): lu:=-(log(b)+4+log(1/sqrt(32)))/(log(a)+4+log(1/sqrt(32))): print(lu): print(`recall that a is `, a1, `and b is `, b1): lu:=evalf(subs({b=b1,a=a1},lu),20): print(`and delta is `): print(``): lprint(lu): print(``): print(`implying the rigorous, but crude upper bound for the irrationality measure, 1+1/delta that equals `): print(``): lprint(evalf(1+1/lu,20)): print(``): print(`confirming the first part of the statement in the title `): print(``): print(`We now need a harder lemma, that we confess that we can't do, but you the reader may be able to `): print(``): print(`Harder Lemma: There exists a rigorously proved constant K1 such that `): print(``): print(`gcd(C1(n),D1(n))=O(exp(K1*n)) `): print(``): print(`Once we find such a constant K1 the delta can be improved to `): print(``): lu:=-(log(b)+4+log(1/sqrt(32))-K1)/(log(a)+4+log(1/sqrt(32))-K1): print(``): print(`By looking at the smallest value of log(gcd(C1(n),D1(n))/n for n between n=`, M/2, `and `, M, `we estimate that K1 can be taken to be`): print(``): lprint(gu[8]): print(``): print(`and the better delta is `): print(``): lu:=subs({a=a1,b=b1,K1=gu[8]},lu): lu:=evalf(lu,20): print(``): lprint(lu): print(``): print(`that implies an irrationality measure, 1+1/delta that equals `): print(``): lprint(evalf(1+1/lu,20)): print(``): print(`-----------------------------------------`): print(``): print(`This ends this paper that took`, time()-t0, `to generate. `): print(``): end: #BestAB(K,M): Explores the empirical delta for the terms between n=M-10 and n=M, in the rational approximations in the generalized Salikhov integral. Using #Int(1/I*(x-4+2*I)^(A*n1)*(x-4-2*I)^(A*n1)*(x-5)^(A*n1)*(x-6+2*I)^(A*n1)*(x-6-2*I)^(A*n1)/(x^(B*n1+1)*(x-10)^(B*n1+1)),x=4-2*I..4+2*I) #for all 1<=A,B<=K with gcd(A,B)=1 and (8,13).Try: #BestAB(10,500); BestAB:=proc(K,M) local gu,A1,B1,A,B,x,n,aluf, si,halev,t0: t0:=time(): gu:=Int(1/I*(x-4+2*I)^(2*A*n)*(x-4-2*I)^(2*A*n)*(x-5)^(2*A*n)*(x-6+2*I)^(2*A*n)*(x-6-2*I)^(2*A*n)/(x^(2*B*n+1)*(x-10)^(2*B*n+1)),x=4-2*I..4+2*I): print(`Trying to Improve Shalikohov's Upper Bound for the Irrationality Measure of Pi`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`In V. Kh. Shalokhov's beautiful article `): print(`"On the Measure of Irrationality of the Number Pi"`): print(`Mathematical Notes,2010, vol.88, no. 4, 563-573 finds a sequence of rational approximations to Pi based on the integral`): print(``): print(subs({A=3,B=5},gu)): print(``): print(`This lead him, after a fairly painful proof to rigorously prove an upper bound of 7.606308 `): print(``): print(`This inspired us to search for variations where the exponents 3 in the numerator, and 5 in the denominator are replaced by other choices`): print(`A,B, relatively prime between 1 and 10 as well as (A,B)= (8,13) in the following generalization`): print(``): print(gu): print(``): print(`Since a rigorous investigation (even by computer) for each case is painful, it makes sense to first find the empirical deltas in `): print(``): print(`abs(Pi-an/bn)<=CONSTANT/bn^(1+delta) `): print(``): print(` for fairly large n. We call this delta the empirical delta, and take, in each case, the smallest for n between`, M-10, `and ` , M): print(``): print(`This will give us an indication of what values of A and B are promising and worth pursuing. `): print(``): print(`Thanks to the recurrences produced, in each case, by the amazing Almkvist-Zeilberger algorithm, it is very easy to crank-out many terms`): print(``): print(`of these integrals, much faster then a direct computation. `): print(``): print(`For the sake of completeness we will list ALL the empirical deltas, even the negative ones, that do not prove irrationality at all.`): print(``): aluf:=[1,1]: si:=min(op(EmpDelAB(1,1,M))): for A1 from 1 to K do for B1 from 1 to K do if igcd(A1,B1)=1 then halev:=min(op(EmpDelAB(A1,B1,M))): print(`(A,B)=`, (A1,B1), `EmpDel=`, halev): if halev>si then si:=halev: aluf:=[A1,B1]: fi: fi: od: od: A1:=8: B1:=13: halev:=min(op(EmpDelAB(8,13,M))): print(`(A,B)=`, (8,13), ` EmpDel= `, halev): if halev>si then si:=halev: aluf:=[A1,B1]: fi: aluf, si: print(`The highest Empirical delta is `, si, `achieved by the case A=`, aluf[1], `and B= `, aluf[2] ): print(`that is better than the empirical delta (A,B)=(3,5) that is`): print(min(op(EmpDelAB(3,5,M)))): print(``): print(`this leads to an irrationality measure around`, evalf(1+1/si,20)): print(``): print(`Hence it is worthwhile to try to investigate this case of`, aluf, `in the hope of beating Salikhov's previous record. `): print(``): print(``): print(`-----------------------------------------`): print(``): print(`This ends this paper that took`, time()-t0, `to generate. `): print(``): end: