print(` Version of April 4, 1997`): print(`This is IRRAT, a Maple package that tries to find`): print(`Apery-Beukers-type irrationality proofs.`): print(`The most current version `): print(`is available at http://www.math.temple.edu/~zeilberg.`): print(`Please report all bugs to: zeilberg@math.temple.edu .`): lprint(``): print(`Written by Doron Zeilberger`): print(`VERY PRELIMINARY!!!!`): ezra:=proc() if args=NULL then print(`This is IRRAT, a Maple packages that tries to find`): print(`Apery-Beukers-type irrationality proofs.`): print(`The most current version `): print(`is available at http://www.math.temple.edu/~zeilberg.`): print(`Please report all bugs to: zeilberg@math.temple.edu .`): lprint(``): print(`Written by Doron Zeilberger`): lprint(``): print(` Contains Procedures `): print(` log2P , QPiP, leg, leg4, Pnm , Pnmtaor, hata7 , Qnm , bestm `): print(` a2rs, Zeta2P , maxb , maxbs , Qm1nm , alpha , delta , bestmP `): print(` delrat, delcon , asymlog2, asymlog20, CATAP, confrec, fApery1 `): print(` For specific help type "ezra(procedure_name)" `): fi: if nops([args])=1 and op(1,[args])=`fApery1` then print(`fApery1(K) finds all recurrences of the form`): print(`n*a(n)+(c1*n+c2)*a(n-1)+(n-1)*a(n-2) that have hope of`): print(`proving irrationality for the limit of the ratio of two solutions`): print(` where -K<=c1,c2<=K`): fi: if nops([args])=1 and op(1,[args])=`confrec` then print(`confrec(ope,n,N,INIa,INIb,L): given a recurrence operator`): print(`ope(n,N) and initial conditions INIa INIb (lists of integers)`): print(`whose size is the same as the Order, and an intger L>10`): print(`outputs a[2*L]/b[2*L] and the corresponding deltas (from delrat)`): print(` for a[L]/b[L], ...., a[L+7]/b[L+7]. `): print(`Here a[n], b[n] are the solutions of the recurrrence `): print(` ope(n,N) x[n]=0 under the initial conditions INIa and INIb, resp.`): fi: if nops([args])=1 and op(1,[args])=`QPiP` then print(`QPiP(P,x) gives the diophantine appx. of Pi/4 and the`): print(` delta corresponding to `): print(`int(int(P/(1+x^2),x=0..1) where P is a polynomial`): print(`in x, with monomials x^m s.t. m==0 mod 2`): fi: if nops([args])=1 and op(1,[args])=`leg2` then print(`leg2(x,n) gives the Legendre-like polynomials, `): print(`D^(2n)[x^(2n)*(1-x^2)^n] over [0,1]`): fi: if nops([args])=1 and op(1,[args])=`leg4` then print(`leg4(x,n) gives the Legendre-like polynomials, `): print(`D^(4n)[x^(4n)*(1-x^4)^n] over [0,1]`): fi: if nops([args])=1 and op(1,[args])=`CATAP` then print(`CATAP(P,x,y) gives the diophantine appx. and the`): print(` delta corresponding to `): print(`int(int(P/(1+x^2*y^2),x=0..1,y=0..1) where P is a polynomial`): print(`in x,y with monomials x^m*y^n s.t. m-n==0 mod 4`): fi: if nops([args])=1 and op(1,[args])=`asymlog20` then print(`asymlog20(K) finds alpha and beta such that`): print(`Pnmtaor(-1,n,0)~alpha^n and int(Pnmtaor(x,n,0)/(1+x),x=0..1)~beta^n`): print(`going empirically up to K. In fact it gives the list of ratios`): print(`and the correspoding delta's s.t. (den(appx))^(1+delta)=1/error`): print(`it is assumed that the gcd of the coeffs. of the polynomials is 1`): print(`(i.e. we don't take advantage of the Hata lemma)`): fi: if nops([args])=1 and op(1,[args])=`asymlog2` then print(`asymlog2(rati,K) finds alpha and beta such that`): print(`Pnmtaor(-1,n,rati*n)~alpha^n and `): print(`int(Pnmtaor(x,n,rati*n)/(1+x),x=0..1)~beta^n`): print(`going empirically up to K. In fact it gives the list of ratios`): print(`and the correspoding delta's s.t. (den(appx))^(1+delta)=1/error`): print(`it is assumed that the gcd of the coeffs. of the polynomials is 1`): print(`(i.e. we don't take advantage of the Hata lemma)`): fi: if nops([args])=1 and op(1,[args])=`delcon` then print(`declcon(alpha), given a constant alpha finds the deltas (as in`): print(`delrat) of the convergents of its continued fraction`): fi: if nops([args])=1 and op(1,[args])=`delrat` then print(`delrat(alpha,rati): Given a constant alpha, and a proposed rational`): print(`appx. rati finds the delta such that`): print(` |alpha-rati|=1/denom(rati)^(1+delta)`): fi: if nops([args])=1 and op(1,[args])=`bestmP` then print(`bestmP(n): Given a positive integer n, finds the m and k such that`): print(`log2P(Pnm(x,n,m),x) is highest, with the corresponding delta`): print(`Here Pnm(x,n,m) is D^n(x^(n-m)*(1-x)^(n+m)))/n!`): fi: if nops([args])=1 and op(1,[args])=`delta` then print(`delta(beta,K): Given a rational number beta, finds the`): print(` constant delta`): print(`in the sequence of appx. to log(2) given by Qnm(x,n,m) with`): print(`n=((1+beta)/2)*n1, m=beta*n1 and n1 ranging over multiples of`): print(`the denominator of beta/2. Recall that the appx. are`): print(`int((1/P(-1)*(P(x)-P(-1))/(1+x), x=0..1) and of |log(2)-p_n/q_n|=`): print(`O(1/q_n^(1+delta))`): fi: if nops([args])=1 and op(1,[args])=`alpha` then print(`alpha(beta,K): Given a rational number beta, finds the constant`): print(` alpha such that Q_{((beta+1)/2)*n,beta*n} ~alpha^(n) up to`): print(` power-law corrections.`): print(`In fact it finds the sequence of appx. up to K`): fi: if nops([args])=1 and op(1,[args])=`Qm1nm` then print(`Qm1nm(n,m) gives D^m(x^n*(1-x)^n)/m! evaluated at x=-1`): fi: if nops([args])=1 and op(1,[args])=`maxbs` then print(`maxbs(bete) finds the maximum of`): print(` f(x):=(x*(1-x))^((1+beta)/2)/(1+x)^beta in the interval [0,1]`): print(` symbolically in beta`): fi: if nops([args])=1 and op(1,[args])=`maxb` then print(`maxb(beta) finds the maximum of`): print(` f(x):=(x*(1-x))^((1+beta)/2)/(1+x)^beta in the interval [0,1]`): fi: if nops([args])=1 and op(1,[args])=`a2rs` then print(`a2rs(r,s) finds the integral int(x^r*y^s/(1-xy),x=0..1,y=0..1)`): print(`Zeta2 is a global constant denoting Zeta(2)`): fi: if nops([args])=1 and op(1,[args])=`bestm` then print(`bestm(n): Given a positive integer n, finds the m such that`): print(`log2P(Qnm(x,n,m),x) is highestm, with the corresponding delta`): print(`Here Qnm(x,n,m) is D^m(x^n*(1-x)^n))/m!`): fi: if nops([args])=1 and op(1,[args])=`Qnm` then print(`Qnm(x,n,m) gives D^m(x^n*(1-x)^n)/m!`): fi: if nops([args])=1 and op(1,[args])=`hata7` then print(`hata7(n): finds the appx. and exponent delta for Hata's appx. to`): print(`ln(2) using P_{7n,n}`): fi: if nops([args])=1 and op(1,[args])=`Pnmtaor` then print(`Pnmtaor(x,n,m): the generalized Legendre polynomials over [0,1]`): fi: if nops([args])=1 and op(1,[args])=`leg` then print(`leg(x,n) gives the Legendre polynomials over [0,1]`): fi: if nops([args])=1 and op(1,[args])=`Pnm` then print(`Pnm(x,n,m) gives the generalized Legendre polynomials over [0,1]`): print(`divided by the gcd of the coefficients`): fi: if nops([args])=1 and op(1,[args])=`log2P` then print(`log2P(P,x): given a polynomial P(x) in x , `): print(`finds int(P/(1+x),x=0..1)/P(-1)`): print(`the error parameter delta`): print(`in the diophantine appx. to log(2), and the rational number`): print(`c:=int((P(x)-P(-1))/(1+x),x=0..1)/P(-1)`): print(` here err=1/denom(c)^(1+delta)`): fi: if nops([args])=1 and op(1,[args])=`Zeta2P` then print(`Zeta2P(P,x,y): given a polynomial P(x,y) in x,y `): print(`finds the diophantine appx.,c, to Zeta(2) implied by performing`): print(`int(P(x,y)/(1-xy),x=0..1,y=0..1)`): print(` followed by the error parameter delta`): print(`in the diophantine appx. to Zeta(2), and the rational number`): print(``): print(` here err=1/denom(c)^(1+delta)`): fi: end: Digits:=50: log2P:=proc(P,x) local c,err,b: c:=int((P-subs(x=-1,P))/(1+x),x=0..1): c:=-c/subs(x=-1,P): err:=evalf(abs(log(2)-c)): b:=denom(c): c,-log(err)/evalf(log(b))-1: end: #delrat(alpha,rati): Given a constant alpha, and a proposed rational #appx. rati finds the delta such that |alpha-rati|=1/denom(rati)^(1+delta) delrat:=proc(alpha,rati) local err,b: if not type(rati,rational) then ERROR(`The second argument to delrat should be a rational number`): fi: err:=evalf(abs(alpha-rati)): b:=denom(rati): -log(err)/evalf(log(b))-1: end: #declcon(alpha), given a constant alpha finds the deltas (as in #delrat) of the convergents of its continued fraction delcon:=proc(alpha) local i,lu,ka: convert(evalf(alpha),confrac,ka): lu:=[]: for i from 3 to nops(ka)/2 do lu:=[op(lu),delrat(alpha,ka[i])]: od: lu: end: #a2rs(r,s) finds the integral int(x^r*y^s/(1-xy),x=0..1,y=0..1) #Zeta2 is a global constant denoting Zeta(2) a2rs:=proc(r1,s1) local gu,k,r,s: if r1=s1 then r:=r1: gu:=Zeta2: for k from 1 to r do gu:=gu-1/k^2: od: RETURN(gu): fi: r:=max(r1,s1): s:=min(r1,s1): gu:=0: for k from s+1 to r do gu:=gu+1/k: od: gu:=gu/(r-s): gu: end: Zeta2P:=proc(P,x,y) local gu,c,err,b,r,s,ar,ars: gu:=0: for r from 0 to degree(P,x) do ar:=coeff(expand(P),x,r): for s from 0 to degree(ar,y) do ars:=coeff(expand(ar),y,s): gu:=expand(gu+ars*a2rs(r,s)): od: od: gu:=expand(gu/coeff(gu,Zeta2,1)): c:=-coeff(gu,Zeta2,0): err:=subs(Zeta2=evalf(Zeta(2)),gu): err:=abs(err): b:=denom(c): c,-log(err)/evalf(log(b))-1: end: #leg(x,n) gives the Legendre polynomials over [0,1] leg:=proc(x,n): x^n*(1-x)^n: diff(",x$n)/n!: expand("): sort("): end: #leg4(x,n) gives the Legendre-like polynomials, #D^(4n)[x^(4n)*(1-x^4)^n] over [0,1] leg4:=proc(x,n): x^(4*n)*(1-x^4)^n: diff(",x$(4*n))/(4*n)!: expand("): sort("): end: #leg2(x,n) gives the Legendre-like polynomials, #D^(2n)[x^(2n)*(1-x^2)^n] over [0,1] leg2:=proc(x,n): x^(2*n)*(1-x^2)^n: diff(",x$(2*n))/(2*n)!: expand("): sort("): end: #Pnm(x,n,m) gives the generalized Legendre polynomials over [0,1] Pnm:=proc(x,n,m) local i,gu,gc: x^(n-m)*(1-x)^(n+m): diff(",x$n)/n!: gu:=expand("): gc:=coeff(gu,x,0): for i from 1 to degree(gu,x) do gc:=gcd(gc,coeff(gu,x,i)): od: gu:=gu/gc: sort(expand(x^m*gu)): end: #P2nm(x,n,m) gives the generalized Legendre leg2, polynomials over [0,1] P2nm:=proc(x,n,m) local i,gu,gc: expand((x)^(2*n-2*m)*(1-x^2)^(n+m)): diff(",x$(2*n))/(2*n)!: gu:=expand("): gc:=coeff(gu,x,0): for i from 1 to degree(gu,x) do gc:=gcd(gc,coeff(gu,x,i)): od: gu:=gu/gc: sort(expand(x^m*gu)): end: #Pnmtaor(x,n,m) gives the generalized Legendre polynomials over [0,1] Pnmtaor:=proc(x,n,m) local i,gu,gc: x^(n-m)*(1-x)^(n+m): diff(",x$n)/n!: gu:=expand("): sort(expand(x^m*gu)): end: #Pnmk(x,n,m,k) gives the generalized Legendre polynomials over [0,1] Pnmk:=proc(x,n,m,k) local i,gu,gc: x^(n-m)*(1-x)^(n+m): diff(",x$n)/n!: gu:=expand(x^k*"): gc:=coeff(gu,x,0): for i from 1 to degree(gu,x) do gc:=gcd(gc,coeff(gu,x,i)): od: gu:=gu/gc: sort(gu): end: #Qnm(x,n,m) gives D^m(x^n*(1-x)^n)/m! with the smallest coeff. 1 Qnm:=proc(x,n,m) local i,gu,gc: gu:=x^(n)*(1-x)^(n): if m>0 then gu:=diff(gu,x$m)/m!: fi: gu:=expand("): gc:=coeff(gu,x,0): for i from 1 to degree(gu,x) do gc:=gcd(gc,coeff(gu,x,i)): od: gu:=gu/gc: sort(gu): end: #Q1nm(x,n,m) gives D^m(x^n*(1-x)^n)/m! Q1nm:=proc(x,n,m) local gu: x^(n)*(1-x)^(n): diff(",x$m)/m!: gu:=expand("): sort(gu): end: #Qm1nmleat(n,m) gives D^m(x^n*(1-x)^n)/m! evaluated at x=-1, slowly Qm1nmleat:=proc(n,m) local x: subs(x=-1,Q1nm(x,n,m)): end: #Qm1nm(n,m) gives D^m(x^n*(1-x)^n)/m! evaluated at x=-1 #apart from sign Qm1nm:=proc(n,m) local gu,k: gu:=0: for k from 0 to n do gu:=gu+binomial(n,k)*binomial(n+k,m): od: gu: end: #hata7(n): finds the appx. and exponent delta for Hata's appx. to #ln(2) using P_{7n,n} hata7:=proc(n) local x: log2P(x^n*Pnm(x,7*n,n),x): end: #bestm(n): Given a positive integer n, finds the m such that #log2P(Qnm(x,n,m),x) is highestm, with the corresponding delta #Here Qnm(x,n,m) is D^m(x^n*(1-x)^n))/m! bestm:=proc(n) local x,m,d1,dm,aluf: d1:=log2P(Qnm(x,n,1),x)[2]: aluf:=1: for m from 2 to trunc((3/2)*n)-1 do dm:=log2P(Qnm(x,n,m),x)[2]: if dm>d1 then d1:=dm: aluf:=m: fi: od: aluf,d1: end: #bestmP(n): Given a positive integer n, finds the m and k such that #log2P(Pnm(x,n,m),x) is highest, with the corresponding delta #Here Pnm(x,n,m) is D^n(x^(n-m)*(1-x)^(n+m)))/n! bestmP:=proc(n) local x,m,d1,dm,aluf: d1:=log2P(Pnm(x,n,1),x)[2]: aluf:=1: for m from 2 to n do dm:=log2P(Pnm(x,n,m),x)[2]: if dm>d1 then d1:=dm: aluf:=m: fi: od: aluf,d1: end: maxb:=proc(beta) local f,x,mu,i,ka,lu: f:=(x*(1-x))^((1+beta)/2)/(1+x)^beta: mu:=diff(log(f),x): mu:=normal(mu): mu:=numer(mu): lu:=[fsolve(mu,x)]: ka:=[]: for i from 1 to nops(lu) do if op(i,lu)>=0 and op(i,lu)<=1 then ka:=[op(ka),abs(subs(x=op(i,lu),f))]: fi: od: max(op(ka)): end: #maxbs(bete) finds the maximum of f(x):=(x*(1-x))^((1+beta)/2)/(1+x)^beta #in the interval [0,1] symbolically in beta maxbs:=proc(beta) local f,x,mu: f:=(x*(1-x))^((1+beta)/2)/(1+x)^beta: mu:=diff(log(f),x): mu:=normal(mu): mu:=numer(mu): subs(x=solve(mu,x)[1],f): end: #alpha(beta,K): Given a rational number beta, finds the constant alpha #such that Q_{((beta+1)/2)*n,beta*n} ~alpha^(n) up to power-law corrections #In fact it finds the sequence of appx. up to K alpha:=proc(beta,K) local q,n1,lu,ku,ku1: if not type(beta,rational) then ERROR(`the first argument must be a rational number`): fi: if not type(K,integer) then ERROR(`the second argument must be an integer`): fi: q:=denom((1+beta)/2): lu:=[]: ku:=Qm1nm(q*(1+beta)/2,beta*q): for n1 from 2 to trunc(K/q) do ku1:=Qm1nm((n1*q)*(1+beta)/2,n1*q*beta): lu:=[op(lu),evalf(ku1/ku)^(1/q)]: ku:=ku1: od: lu,op(nops(lu),lu): end: #delta(beta,K): Given a rational number beta, finds the constant delta #in the sequence of appx. to log(2) given by Qnm(x,n,m) with #n=((1+beta)/2)*n1, m=beta*n1 and n1 ranging over multiples of #the denominator of beta/2. Recall that the appx. are #int((1/P(-1)*(P(x)-P(-1))/(1+x), x=0..1) and of |log(2)-p_n/q_n|= #O(1/q_n^(1+delta)) delta:=proc(beta,K) local beta1,alpha1: beta1:=1/maxb(beta): alpha1:=alpha(beta,K)[2]: (log(alpha1)+log(beta1))/(1.+log(alpha1)): evalf("-1): end: #asymlog2(rati,K) finds alpha and beta such that #Pnmtaor(-1,n,rati*n)~alpha^n and #int(Pnmtaor(x,n,rati*n)/(1+x),x=0..1)~beta^n #going empirically up to K. In fact it gives the list of ratios #and the correspoding delta's such that (den(appx))^(1+delta)=1/error #it is assumed that the gcd of the coeffs. of the polynomials is 1 #(i.e. we don't take advantage of the Hata lemma) asymlog2:=proc(rati,K) local alpha,beta,alpha1,beta1,delta,delta1,q,gu1,gu,lu1,lu,P,x: if not type(rati,rational) then ERROR(`The argument of asymlog2 must be a rational number`): fi: q:=denom(rati): alpha:=[]: beta:=[]: delta:=[]: P:=Pnmtaor(x,q,rati*q): gu:=subs(x=-1,P): lu:=abs(evalf(int(P/(1+x),x=0..1))): for n from 2*q by q to K do P:=Pnmtaor(x,n,rati*n): gu1:=subs(x=-1,P): lu1:=abs(evalf(int(P/(1+x),x=0..1))): alpha1:=evalf((gu1/gu)^(1/q)): beta1:=evalf((lu1/lu)^(1/q)): delta1:=(log(alpha1)-log(beta1))/(log(alpha1)+1+rati)-1: alpha:=[op(alpha),alpha1]: beta:=[op(beta),beta1]: delta:=[op(delta),delta1]: gu:=gu1: lu:=lu1: od: alpha,beta,delta: end: #asymlog20(K) finds alpha and beta such that #Pnmtaor(-1,n,0)~alpha^n and int(Pnmtaor(x,n,0)/(1+x),x=0..1)~beta^n #going empirically up to K. In fact it gives the list of ratios #and the correspoding delta's such that (den(appx))^(1+delta)=1/error #it is assumed that the gcd of the coeffs. of the polynomials is 1 #(i.e. we don't take advantage of the Hata lemma) asymlog20:=proc(K) local alpha,beta,alpha1,beta1,delta, delta1, q,gu1,gu,lu1,lu,P,x: delta:=[]: alpha:=[]: beta:=[]: P:=Pnmtaor(x,1,0): gu:=subs(x=-1,P): lu:=abs(evalf(int(P/(1+x),x=0..1))): for n from 2 to K do P:=Pnmtaor(x,n,0): gu1:=subs(x=-1,P): lu1:=abs(evalf(int(P/(1+x),x=0..1))): alpha1:=evalf(gu1/gu): beta1:=evalf(lu1/lu): delta1:=(log(alpha1)-log(beta1))/(log(alpha1)+1)-1: alpha:=[op(alpha),alpha1]: beta:=[op(beta),beta1]: delta:=[op(delta),delta1]: gu:=gu1: lu:=lu1: od: alpha,beta,delta: end: #C2rs(r,s) finds the integral int(x^r*y^s/(1+x^2*y^2),x=0..1,y=0..1) #in terms of a combination of Catalan's constant, to be denoted by #the global variable CATA, and rational numbers #r and s should be non-negative integers s.t. r-s==0 mod 4. If this #does not hold, and error message ensues C2rs:=proc(r1,s1) local gu,k,r,s,R,m: if not type(r1/2,integer) or not type(s1/2,integer) or r1<0 or s1<0 or r1-s1 mod 4<>0 then ERROR(`The inputs should be non-neg. even integ. whose diff. is div. by 4`): fi: if r1=s1 then if type(r1/4,integer) then m:=r1/4: gu:=CATA: for k from 0 to 2*m-1 do gu:=gu-(-1)^k/(2*k+1)^2: od: RETURN(gu): fi: if type((r1-2)/4,integer) then m:=(r1-2)/4: gu:=CATA: for k from 0 to 2*m do gu:=gu-(-1)^k/(2*k+1)^2: od: RETURN(-gu): fi: fi: r:=max(r1,s1): s:=min(r1,s1): gu:=0: R:=(r-s)/4: for k from 0 to 2*R-1 do gu:=gu+(-1)^k/(2*k+s+1): od: gu:=gu/4/R: gu: end: #CATAP(P,x,y) gives the diophantine appx. and the # delta corresponding to #int(int(P/(1+x^2*y^2),x=0..1,y=0..1) where P is a polynomial #in x,y with monomials x^m*y^n s.t. m-n=mod 4 CATAP:=proc(P,x,y) local gu,c,err,b,r,s,ar,ars: gu:=0: for r from 0 to degree(P,x) do ar:=coeff(expand(P),x,r): for s from 0 to degree(ar,y) do ars:=coeff(expand(ar),y,s): if not (r-s mod 4=0) and ars<>0 then print(r-s mod 4): ERROR(`The polynomial P contains the monomial`, x^r*y^s, `which is no good`): fi: if (r=1 mod 2 or s=1 mod 2) and ars<>0 then ERROR(`The polynomial P contains the monomial`, x^r*y^s, `which is no good`): fi: if ars<>0 then gu:=expand(gu+ars*C2rs(r,s)): fi: od: od: if coeff(gu,CATA,1)=0 then print(`The polynomial`,P,`does not yield an appx.`): RETURN(0): fi: gu:=expand(gu/coeff(gu,CATA,1)): c:=-coeff(gu,CATA,0): err:=subs(CATA=evalf(Catalan),gu): err:=abs(err): b:=denom(c): if b<>1 then RETURN(c,-log(err)/evalf(log(b))-1): else RETURN(c): fi: end: #kak(r,s) checks the validity of procedure C2rs kak:=proc(r,s) local gu,x,y: gu:=C2rs(r,s): gu:=subs(CATA=evalf(Catalan),gu): gu:=gu-evalf(Int(Int(x^r*y^s/(1+x^2*y^2),x=0..1),y=0..1)): gu: end: #QPir(r) finds the integral int(x^r/(1+x^2),x=0..1,y=0..1) #in terms of a combination of Pi/4, to be denoted by #the global variable QPi, and rational number. # #r should be non-negative integers s.t. r==0 mod 2. If this #does not hold, and error message ensues QPir:=proc(r1) local gu,k,r,s,R,m: if not type(r1/2,integer) or r1<0 then ERROR(`The inputs should be non-neg. even integ.`): fi: if type(r1/4,integer) then m:=r1/4: gu:=QPi: for k from 0 to 2*m-1 do gu:=gu-(-1)^k/(2*k+1): od: RETURN(gu): fi: if type((r1-2)/4,integer) then m:=(r1-2)/4: gu:=QPi: for k from 0 to 2*m do gu:=gu-(-1)^k/(2*k+1): od: RETURN(-gu): fi: end: #QPi(P,x) gives the diophantine appx. of Pi/4 and the # delta corresponding to #int(int(P/(1+x^2),x=0..1) where P is a polynomial #in x, with monomials x^m s.t. m==0 mod 2 QPiPold:=proc(P,x) local gu,c,err,b,r,s,ar,ars: gu:=0: for r from 0 to degree(P,x) do ar:=coeff(expand(P),x,r): if not type(r/2,integer) and ar<>0 then ERROR(`The polynomial P contains the monomial`, x^r, `which is no good`): fi: if ar<>0 then gu:=expand(gu+ar*QPir(r)): fi: od: if coeff(gu,QPi,1)=0 then print(`The polynomial`,P,`does not yield an appx.`): RETURN(0): fi: gu:=expand(gu/coeff(gu,QPi,1)): c:=-coeff(gu,QPi,0): err:=subs(QPi=evalf(Pi/4),gu): err:=abs(err): b:=denom(c): if b<>1 then RETURN(c,-log(err)/evalf(log(b))-1): else RETURN(c): fi: end: #QPiP(P,x): More General and simpler version of QPiPold QPiP:=proc(P,x) local gu,c,err,b,r,s,ar,ars,sher,quot: gu:=0: quot:=quo(P,1+x^2,x,sher): if degree(sher,x)=1 then ERROR(`Bad P`): fi: gu:=int(quot,x=0..1)+sher*QPi: if coeff(gu,QPi,1)=0 then print(`The polynomial`,P,`does not yield an appx.`): RETURN(0): fi: gu:=expand(gu/coeff(gu,QPi,1)): c:=-coeff(gu,QPi,0): err:=subs(QPi=evalf(Pi/4),gu): err:=abs(err): b:=denom(c): if b<>1 then RETURN(c,-log(err)/evalf(log(b))-1): else RETURN(c): fi: end: confrec:=proc(ope,n,N,INIa,INIb,L) local mu,ope1,ORDER,i,n1,gadol,a,b,lu,c: if L<10 then ERROR(`L>10`): fi: gadol:=degree(ope,N): ope1:=subs(n=n-gadol,ope)/N^gadol: ope1:=expand(ope1): ORDER:=-ldegree(ope1,N): if nops(INIa)<>ORDER or nops(INIb)<>ORDER then ERROR(`The 4th or 5th arguments should both hve`, ORDER, ` terms `): fi: for i from 0 to ORDER-1 do a[i]:=INIa[i+1]: b[i]:=INIb[i+1]: od: for n1 from ORDER to 3*ORDER do lu:=0: for i from 1 to ORDER do lu:=lu+subs(n=n1,coeff(ope1,N,-i))*a[n1-i]: od: lu:=-lu/subs(n=n1,coeff(ope1,N,0)): a[n1]:=lu: lu:=0: for i from 1 to ORDER do lu:=lu+subs(n=n1,coeff(ope1,N,-i))*b[n1-i]: od: lu:=-lu/subs(n=n1,coeff(ope1,N,0)): b[n1]:=lu: od: for n1 from 3*ORDER+1 to 2*L do lu:=0: for i from 1 to ORDER do lu:=lu+subs(n=n1,coeff(ope1,N,-i))*a[n1-i]: od: lu:=-lu/subs(n=n1,coeff(ope1,N,0)): a[n1]:=lu: lu:=0: for i from 1 to ORDER do lu:=lu+subs(n=n1,coeff(ope1,N,-i))*b[n1-i]: od: lu:=-lu/subs(n=n1,coeff(ope1,N,0)): b[n1]:=lu: od: if b[2*L]<>0 then c:=a[2*L]/b[2*L]: else RETURN(LOTOV): fi: mu:=[]: for i from L to L+7 do mu:=[op(mu),delrat(c,a[i]/b[i])]: od: c,mu: end: #fApery1(K) finds all recurrences of the form #n*a(n)+(c1*n+c2)*a(n-1)+(n-1)*a(n-2) that have hope of #proving irrationality for the limit of the ratio of two solutions # where -K<=c1,c2<=K fApery1:=proc(K) local ope,gu,c1,c2,lu: lu:={}: INIa:=[1,0]: INIb:=[0,1]: for c1 from -K to K do for c2 from -K to K do ope:=n+(c1*n+c2)*N^(-1)+(n-1)*N^(-2): gu:=confrec(ope,n,N,INIa,INIb,30): if gu<>LOTOV then gu1:=gu[1]: gu2:=gu[2]: gu2:=min(op(gu2)): if gu2> 0.05 then lu:=lu union {[ope,gu1,gu2]}: fi: fi: od: od: lu: end: #fApery1t(K) Terse version of fApery1 finds all recurrences of the form #n*a(n)+(c1*n+c2)*a(n-1)+(n-1)*a(n-2) that have hope of #proving irrationality for the limit of the ratio of two solutions # where -K<=c1,c2<=K fApery1t:=proc(K) local ope,gu,c1,c2,lu: lu:={}: INIa:=[1,0]: INIb:=[0,1]: for c1 from -K to K do for c2 from -K to K do ope:=n+(c1*n+c2)*N^(-1)+(n-1)*N^(-2): gu:=confrec(ope,n,N,INIa,INIb,30): if gu<>LOTOV then gu1:=gu[1]: gu2:=gu[2]: gu2:=min(op(gu2)): if gu2> 0.05 then lu:=lu union {ope}: fi: fi: od: od: lu: end: #fApery1ts(K,c1) Terse version of fApery1 finds all recurrences of the form #n*a(n)+(c1*n+c2)*a(n-1)+(n-1)*a(n-2) that have hope of #proving irrationality for the limit of the ratio of two solutions # where -K<=c2<=K fApery1ts:=proc(K,c1) local ope,gu,c2,lu: lu:=[]: INIa:=[1,0]: INIb:=[0,1]: for c2 from -K to K do ope:=n+(c1*n+c2)*N^(-1)+(n-1)*N^(-2): gu:=confrec(ope,n,N,INIa,INIb,30): if gu<>LOTOV then gu1:=gu[1]: gu2:=gu[2]: gu2:=min(op(gu2)): if gu2> 0.05 then lu:=[op(lu),ope]: fi: fi: od: lu: end: fApery1T:=proc(K) local lu,c1: lu:=[]: for c1 from -K to K do lu:=[op(lu),fApery1ts(K,c1)]: od: lu: end: