###################################################################### ##AperyAppx: Save this file as AperyAppx. To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read AperyAppx : # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg@math.rutgers.edu. # ###################################################################### #Created: April 22, 2002 #This version: April 22, 2002 #AperyAppx: A Maple package to study Rational and #approximations inspired by Apery's Original Method #As described in:"Interpolation de Fractions Continues et # Irrationalite de Certaines Constantes": #Bull. de la section des sci. du C.T.H.S. no. 3, p. 37-53 #Bib. Nat. Paris, Math Reviews 83b:10039 #It is one of the packages accompanying Doron Zeilberger's #article COMPUTERIZED DECONSTRUCTION, available from #his website #Please report bugs to zeilberg@math.rutgers.edu print(` AperyAppx: A Maple package to study Rational Approximations`): print(` inspired by Apery's original method`): print(`As described in his masterpiece:`): print(`"Interpolation de Fractions Continues et Irrationalite de `): print(` Certaines Constantes":`): print(`Bull. de la section des sci. du C.T.H.S. no. 3, p. 37-53`): print(` (Math Reviews 83b:10039) .`): print( ` It is one of the packages accompanying Doron Zeilberger's`): print(` article COMPUTERIZED DECONSTRUCTION, available from `): print(` his website `): print(``): print(`The most current version of this program, and article, `): print(` are available from Zeilberger's homepage`): print(` http://www.math.rutgers.edu/~zeilberg/ ,`): print(` by clicking at the appropriate places. `): print(``): print(` Please report bugs to zeilberg@math.rutgers.edu `): print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`): print(`For a list of the MAIN procedures, type: ezra(); `): print(``): print(`For a list of ALL procedures, type: ezra1(); `): print(`For help with a specific procedrue, type: ezra(procedure_name); `): ezra1:=proc(): if args=NULL then print(` All procedures are: AperyB, delt, `): print(` DeltaSeq, DeltaSeq1, FindOpe, FindOpeNew, NumerAppx,`): print(` OneStep, OneStep1, Pol, RatAppx,RatAppx1, `): print(` RatAppxOld, RatAppxOld1, RatAppxQ1, RatAppxQ `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` The Main procedures are: AperyB, delt, DeltaSeq1, NumerAppx, `): print(` OneStep, RatAppx `): print(``): print(`With help with a procdure, type: ezra(procedure_name);`): print(`For example, for help with AperyB, type: ezra(AperyB); `): fi: if nops([args])=1 and op(1,[args])=AperyB then print(`AperyB(pol,n,h,degh): Given a polynomial pol(n) tries to find`): print(`the accelerating B(n,h) of degree degh in h`): print(`such that P_h/Q_h->`): print(`(B(n,h)*P(n)+pol(n+1)*P(n+1)+Q(n+1))/(B(n,h)*Q(n)+pol(n+1)*Q(n+1))`): print(`is an improvement of sum(1/pol(i),i=n+1..infinity)`): print(`Returns 0 if it fails`): print(`For example, try: AperyB(n^2,n,h,2); AperyB(n^3,n,h,3);`): fi: if nops([args])=1 and op(1,[args])=delt then print(`delt(a,c): Given a rational number a and a constant`): print(`c, finds the delta such that |a-c|=1/denom(a)^(1+delta)`): print(`For example, try delta(22/7,evalf(Pi)):`): fi: if nops([args])=1 and op(1,[args])=DeltaSeq then print(`DeltaSeq(pol,Cf,n,L): The sequence of deltas for the diophantine`): print(`approximations offered by NumerAppx to`): print(`sum(Cf(i)/pol(i),i=1..infinity)`): fi: if nops([args])=1 and op(1,[args])=DeltaSeq1 then print(`DeltaSeq1(pol,Cf,n,L,Const): `): print(`The sequence of deltas for the diophantine`): print(`approximations offered by NumerAppx to`): print(`Const=sum(Cf(i)/pol(i),i=1..infinity)`): print(`For example, try: DeltaSeq1(-(2*n-1)^2,(-1)^n,n,20,Catalan); `): fi: if nops([args])=1 and op(1,[args])=FindOpe then print(`FindOpe(Q,d,n,N): tries to find a second-order operator of the`): print(`form ope(N,n)=P2*N^2+P1(n)*N+P0(n) where P0,P1,P2 are`): print(`polynomials of degree d`): print(`annihilating the polynomial Q(n)`): print(`For example, try FindOpe(n,0,n,N);`): fi: if nops([args])=1 and op(1,[args])=FindOpeNew then print(`FindOpeNew(Q,d,n,N,pol): tries to find a second-order`): print(` operator of the`): print(`form ope(N,n)=pol(n+2)*P2*N^2+P1(n)*N+pol(n+1)*P0(n) `): print(`where ope(N,n) is`): print(`of degree d in n,`): print(`annihilating the polynomial Q(n)`): fi: if nops([args])=1 and op(1,[args])=NumerAppx then print(`NumerAppx(pol,Cf,n,d2,L): The numerical appx`): print(` for sum(Cf(n)/pol(n),n=1..infinity) offered by`): print(`#sum(Cf(n)/pol(n),n=1..L) +R_d2(L)*Cf(L)`): print(`where R_d2(n) is RatAppx(pol,Cf,n,d2);`): print(`For example, try: NumerAppx(1,1/n!^2,n,5,5);`): fi: if nops([args])=1 and op(1,[args])=OneStep then print(`OneStep1(R,pol,n,L): Given a rational function R(n)(=P(n)/Q(n)) `): print(`such that`): print(`R(n)-R(n-1)+1/pol(n) is small, finds a better one such that`): print(`R1(n):=`): print(`(B(n)*P(n)+pol(n+1)*P(n+1)+Q(n+1))/(B(n)*Q(n)+pol(n+1)*Q(n+1))`): print(`the degree of the numerator of R1(n)-R1(n-1)+1/pol(n) is as small`): print(`as possible(where the degree of B (in n) is that of pol(n))`): print(`The output is the good B followed by the modified R`): print(`For example, try: OneStep(0,n^3,n);`): fi: if nops([args])=1 and op(1,[args])=OneStep1 then print(`OneStep1(R,pol,n,L): Given a rational function R(n)(=P(n)/Q(n)) `): print(`such that`): print(`R(n)-R(n-1)+1/pol(n) is small, finds a better one such that`): print(`R1(n):=`): print(`(B(n)*P(n)+pol(n+1)*P(n+1)+Q(n+1))/(B(n)*Q(n)+pol(n+1)*Q(n+1))`): print(`the degree of the numerator of R1(n)-R1(n-1)+1/pol(n) is L`): print(`(where the degree of B (in n) is that of pol(n))`): print(`For example, try: OneStep1(0,n^3,n,0);`): fi: if nops([args])=1 and op(1,[args])=Pol then print(`Pol(List,h,h0,deg): Given a list of expressions, List,`): print(`guesses a polynomial P,of degree deg such that`): print(`it fits List=[P(h0),P(h0+1),P(h0+nops(List)-1)], and returns`): print(`0 if none exists `): fi: if nops([args])=1 and op(1,[args])=QuasiRatAppxOld1 then print(`QuasiRatAppxOld1(pol,n,deg,m,L): Tries to find the `): print(`Quasi-rational function of period m`): print(`R(n), with monic denominator of of degree deg,`): print(`R(n)-R(n-1)+1/pol(n) has numerator of degree L`): print(`For example, try: QuasiRatAppxOld1(n^2,n,3,1,0);`): fi: if nops([args])=1 and op(1,[args])=RatAppx then print(`RatAppx(pol,Cf,n,deg): Tries to find the rational function`): print(` R(n), with monic denominator of of degree deg, `): print(` (Cf(n)*R(n)-Cf(n-1)*R(n-1)+Cf(n)/pol(n))/Cf(n)`): print(` has numerator of degree as small as possible `): print(` For example, try: RatAppx(n^2,1,n,3); `): fi: if nops([args])=1 and op(1,[args])=RatAppx1 then print(`RatAppx1(pol,Cf,n,deg,L): Tries to find the rational function`): print(` R(n), with monic denominator of of degree deg, `): print(` (Cf(n)*R(n)-Cf(n-1)*R(n-1)+Cf(n)/pol(n))/Cf(n)`): print(` has numerator of degree L `): print(` For example, try: RatAppx1(n^2,1,n,3,0); `): fi: if nops([args])=1 and op(1,[args])=RatAppxOld then print(`RatAppxOld(pol,n,deg): Tries to find the rational function`): print(` R(n), with monic denominator of of degree deg, `): print(` R(n)-R(n-1)+1/pol(n) has numerator of degree `): print(` as small as possible `): print(` For example, try: RatAppxOld(n^2,n,3); `): fi: if nops([args])=1 and op(1,[args])=RatAppxOld1 then print(`RatAppxOld1(pol,n,deg,L): Tries to find the rational function`): print(`R(n), with monic denominator of of degree deg,`): print(` R(n)-R(n-1)+1/pol(n) has numerator of degree L `): print(` For example, try: RatAppxOld1(n^2,n,3,0);`): fi: end: #RatAppxOld1(pol,n,deg,L): Tries to find the rational function #R(n), with monic denominator of of degree deg, #R(n)-R(n-1)+1/pol(n) has numerator of degree L #For example, try: RatAppxOld1(n^2,n,3,0); RatAppxOld1:=proc(pol,n,d2,L) local mone,mekh,a,b,eq,var,i,gu,R,d1,Deg1: Deg1:=degree(numer(pol),n)-degree(denom(pol),n): d1:=d2-(Deg1-1): var:={}: mone:=0: for i from 0 to d1 do mone:=mone+a[i]*n^i: var:=var union {a[i]}: od: mekh:=n^d2: for i from 0 to d2-1 do mekh:=mekh+b[i]*n^i: var:=var union {b[i]}: od: R:=mone/mekh: gu:=expand(numer(normal(R-subs(n=n-1,R)+1/pol))): eq:={}: for i from L+1 to degree(gu,n) do eq:=eq union {coeff(gu,n,i)}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: if nops({var})>1 then print(`Too many solutions`): RETURN(0): fi: normal(subs(var,R)): end: #RatAppxOld(pol,n,deg): Tries to find the rational function #R(n), with monic denominator of of degree deg, #R(n)-R(n-1)-1/pol(n) has numerator of degree #as small as possible #For example, try: RatAppxOld(n^2,n,3); RatAppxOld:=proc(pol,n,d2) local L,gu: if d2=degree(pol,n)-1`): fi: gu:=0: for L from 0 while gu=0 do gu:=RatAppxOld1(pol,n,d2,L): if gu<>0 then RETURN(gu): fi: od: end: #OneStep1(R,pol,n,L): Given a rational function R(n)(=P(n)/Q(n)) #such that #R(n)-R(n-1)+1/pol(n) is small, finds a better one such that #R1(n):= #(B(n)*P(n)+pol(n+1)*P(n+1)+Q(n+1))/(B(n)*Q(n)+pol(n+1)*Q(n+1)) #the degree of the numerator of R1(n)-R1(n-1)+1/pol(n) is L #(where the degree of B (in n) is that of pol(n)) #For example, try: OneStep1(0,n^3,n,0); OneStep1:=proc(R,pol,n,L) local P,Q,b,B,eq,var,i,R1,lu,var1,ulai,R1a,Ba,shi,alufB,alufR: _EnvAllSolutions:=true: P:=numer(R): Q:=denom(R): B:=0: var:={}: for i from 0 to degree(pol,n) do B:=B+b[i]*n^i: var:=var union {b[i]}: od: R1:=(B*P+subs(n=n+1,pol*P+Q))/(B*Q+subs(n=n+1,pol*Q)): lu:=expand(numer(normal(R1-subs(n=n-1,R1)+1/pol))): eq:={}: for i from L+1 to degree(lu,n) do eq:=eq union {coeff(lu,n,i)}: od: var1:=[solve(eq,var)]: if var1=[] then RETURN(0): fi: alufB:=subs(var1[1],B): alufR:=normal(subs(var1[1],R1)): shi:=degree(denom(alufR),n): for i from 2 to nops(var1) do Ba:=subs(var1[i],B): R1a:=normal(subs(var1[i],R1)): ulai:=degree(denom(R1a),n): if ulai>shi then shi:=ulai: alufB:=Ba: alufR:=R1a: fi: od: if degree(denom(alufR),n)=0 then RETURN(0): fi: alufB,alufR: end: #OneStep(R,pol,n): Given a rational function R(n)(=P(n)/Q(n)) #such that #R(n)-R(n-1)+1/pol(n) is small, finds a better one such that #R1(n):= #(B(n)*P(n)+pol(n+1)*P(n+1)+Q(n+1))/(B(n)*Q(n)+pol(n+1)*Q(n+1)) #the degree of the numerator of R1(n)-R1(n-1)+1/pol(n) is as small #as possible (where the degree of B (in n) is that of pol(n)) #The output is the good B followed by the imporved R #For example, try: OneStep(0,n^3,n); OneStep:=proc(R,pol,n) local L,gu: gu:=0: #for L from 0 to max(degree(pol,n)-1,degree(denom(R),n)) while gu=0 do for L from 0 while gu=0 do gu:=OneStep1(R,pol,n,L): if gu<>0 then RETURN(gu): fi: od: 0: end: #Pol(List,h,h0,deg): Given a list of expressions, List, #guesses a polynomial P,of degree deg such that #it fits List=[P(h0),P(h0+1),P(h0+nops(List)-1)], and returns #0 if none exists Pol:=proc(List,h,h0,deg) local pol,eq,var,i,a: var:={}: pol:=0: for i from 0 to deg do pol:=pol+a[i]*h^i: var:=var union {a[i]}: od: eq:={}: for i from 1 to nops(List) do eq:=eq union {expand(subs(h=h0+i-1,pol)-List[i])}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: expand(subs(var,pol)): end: #delt(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)): delt:=proc(a,c): Digits:=50: evalf(-log(abs(c-a))/log(denom(a)))-1: end: #AperyB(pol,n): Given a polynomial pol(n) tries to find #the accelerating B(n,h) of degree degh in h #such that P_h/Q_h-> #(B(n,h)*P(n)+pol(n+1)*P(n+1)+Q(n+1))/(B(n,h)*Q(n)+pol(n+1)*Q(n+1)) #is an improvement of sum(1/pol(i),i=n+1..infinity) #Returns 0 if it fails #For example, try: AperyB(n^2,n,h,2); AperyB(n^3,n,h,3); AperyB:=proc(pol,n,h,degh) local gu,i,lu: gu:=[OneStep(0,pol,n)]: if gu=[0] then RETURN(0): fi: lu:=[]: for i from 1 to degh+2 do lu:=[op(lu),gu]: gu:=[OneStep(gu[2],pol,n)]: if gu=[0] then RETURN(0): fi: od: lu:=[seq(lu[i][1],i=1..nops(lu))]: Pol(lu,h,0,degh): end: #RatAppx1(pol,Cf,n,deg,L): #Given a rational function pol(n), and a closed form sequence #Cf(n), tries to find the rational function #R(n), with monic denominator of of degree deg, #(Cf(n)*R(n)-Cf(n-1)R(n-1)-Cf(n)/pol(n))/Cf(n) has numerator of #degree L; For example, try: RatAppx1(n^2,1,n,3,0); RatAppx1:=proc(pol,Cf,n,d2,L) local mone,mekh,a,b,eq,var,i,gu,R,d1,Deg1,ku: Deg1:=degree(numer(pol),n)-degree(denom(pol),n): d1:=d2-(Deg1-1): var:={}: mone:=0: for i from 0 to d1 do mone:=mone+a[i]*n^i: var:=var union {a[i]}: od: mekh:=n^d2: for i from 0 to d2-1 do mekh:=mekh+b[i]*n^i: var:=var union {b[i]}: od: R:=mone/mekh: ku:=normal(expand(normal(simplify(subs(n=n-1,Cf)/Cf)))): gu:=expand(numer(normal(R-ku*subs(n=n-1,R)+1/pol))): eq:={}: for i from L+1 to degree(gu,n) do eq:=eq union {coeff(gu,n,i)}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: if nops({var})>1 then print(`Too many solutions`): RETURN(0): fi: normal(subs(var,R)): end: #RatAppx(pol,Cf,n,deg): Given a rational function pol(n) #and a closed-form function Cf(n) #Tries to find the rational function #R(n), with monic denominator of of degree deg, #(Cf(n)*R(n)-Cf(n-1)*R(n-1)-Cf(n)/pol(n))/Cf(n) #has numerator of degree #as small as possible #For example, try: RatAppx(n^2,1,n,3); RatAppx:=proc(pol,Cf,n,d2) local L,gu: option remember: if d2=degree(pol,n)-1`): fi: gu:=0: for L from 0 while gu=0 do gu:=RatAppx1(pol,Cf,n,d2,L): if gu<>0 then RETURN(gu): fi: od: end: #NumerAppx(pol,Cf,n,d2,L): The numerical appx # for sum(Cf(n)/pol(n),n=1..infinity) offered by #sum(Cf(n)/pol(n),n=1..L) +R_d2(L)*Cf(L) #where R_d2(n) is RatAppx(pol,Cf,n,d2); #For example, try: NumerAppx(1,1/n!^2,n,5,5); NumerAppx:=proc(pol,Cf,n,d2,L) local R,i: R:=RatAppx(pol,Cf,n,d2): convert([seq(subs(n=i,Cf/pol),i=1..L)],`+`)+subs(n=L,R*Cf): end: #RatAppxQ1(pol,n,Q,L): Tries to find the rational function #R(n)=P/Q, with monic denominator of of degree deg, #R(n)-R(n-1)+1/pol(n) has numerator of degree L #For example, try: RatAppxQ1(n^2,n,n^4,3,0); RatAppxQ1:=proc(pol,n,Q,L) local mone,a,eq,var,i,gu,R,d1,Deg1,d2: d2:=degree(Q,n): Deg1:=degree(numer(pol),n)-degree(denom(pol),n): d1:=d2-(Deg1-1): var:={}: mone:=0: for i from 0 to d1 do mone:=mone+a[i]*n^i: var:=var union {a[i]}: od: R:=mone/Q: gu:=expand(numer(normal(R-subs(n=n-1,R)+1/pol))): eq:={}: for i from L+1 to degree(gu,n) do eq:=eq union {coeff(gu,n,i)}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: if nops({var})>1 then print(`Too many solutions`): RETURN(0): fi: normal(subs(var,R)): end: #RatAppxQ(pol,n,Q): Tries to find the rational function #R(n)=P/Q(n), #R(n)-R(n-1)-1/pol(n) has numerator of degree #as small as possible #For example, try: RatAppxOld(n^2,n,3); RatAppxQ:=proc(pol,n,Q) local L,gu,d2: d2:=degree(Q,n): if d2=degree(pol,n)-1`): fi: gu:=0: for L from 0 while gu=0 do gu:=RatAppxQ1(pol,n,Q,L): if gu<>0 then RETURN(gu): fi: od: end: #NumerAppxQ(pol,n,Q,L): The numerical appx # for sum(1/pol(n),n=1..infinity) offered by #sum(1/pol(n),n=1..L) +R_d2(L) #where R_d2(n) is RatAppxQ(pol,n,Q); #For example, try: NumerAppxQ(n^2,n,n^3,5,5); NumerAppxQ:=proc(pol,n,Q,L) local R,i: R:=RatAppxQ(pol,n,Q): convert([seq(subs(n=i,1/pol),i=1..L)],`+`)+subs(n=L,R): end: #FindOpe(Q,d,n,N): tries to find a second-order operator of the #form ope(N,n)=P2*N^2+P1(n)*N+P0(n) where P0,P1,P2 are #polynomials of degree d #annihilating the polynomial Q(n) #For example, try FindOpe(n,0,n,N); FindOpe:=proc(Q,d,n,N) local ope,P0,P1,P2,eq,var,i,a,b,c,gu: var:={}: P0:=0:P1:=0:P2:=0: for i from 0 to d do P0:=P0+a[i]*n^i: P1:=P1+b[i]*n^i: P2:=P2+c[i]*n^i: var:=var union {a[i],b[i],c[i]}: od: gu:=P2*subs(n=n+2,Q)+P1*subs(n=n+1,Q)+P0*Q; ope:=P2*N^2+P1*N+P0; gu:=expand(gu): eq:={}: for i from 0 to degree(gu,n) do eq:=eq union {coeff(gu,n,i)}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: ope:=subs(var,ope): if ope<>0 then RETURN(yafe(ope,N)): else RETURN(0): fi: end: #FindOpeNew(Q,d,n,N,pol): tries to find a second-order operator of the #form ope(N,n)=pol(n+2)*P2*N^2+P1(n)*N+pol(n+1)*P0(n) where P0,P1,P2 are #polynomials of degree d #annihilating the polynomial Q(n) #For example, try FindOpeNew(n,0,n,N,n); FindOpeNew:=proc(Q,d,n,N,pol) local ope,P0,P1,P2,eq,var,i,a,b,c,gu,i1,j1,ope1: #option remember: var:={}: P0:=0:P1:=0:P2:=0: for i from 0 to d-degree(pol,n) do P0:=P0+a[i]*n^i: var:=var union {a[i]}: od: for i from 0 to d do P1:=P1+b[i]*n^i: var:=var union {b[i]}: od: for i from 0 to d-degree(pol,n) do P2:=P2+c[i]*n^i: var:=var union {c[i]}: od: gu:=P2*subs(n=n+2,Q*pol)+P1*subs(n=n+1,Q)+P0*Q*subs(n=n+1,pol); ope:=P2*subs(n=n+2,pol)*N^2+P1*N+P0*subs(n=n+1,pol); gu:=expand(gu): eq:={}: for i from 0 to degree(gu,n) do eq:=eq union {coeff(gu,n,i)}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: ope:=expand(subs(var,ope)): if ope<>0 then ope1:=yafe(ope,N): ope:=expand(ope1): for i1 from 0 to degree(ope,n) do for j1 from 0 to degree(coeff(ope,n,i1),N) do if not type(coeff(coeff(ope,n,i1),N,j1),numeric) then RETURN(0): fi: od: od: RETURN(ope1): else RETURN(0): fi: end: #yafe(ope,N): given a recurrence operator in the shift #N, outputs it in nice form yafe:=proc(ope,N) local ope1,ope2,i: ope1:=expand(ope): ope1:=ope1/coeff(ope1,N,degree(ope1,N)): ope1:=normal(ope1): ope1:=numer(ope1): ope2:=0: for i from ldegree(ope1,N) to degree(ope1,N) do ope2:=ope2+factor(coeff(ope1,N,i))*N^i: od: ope2: end: #DeltaSeq(pol,Cf,n,L): The sequence of deltas for the diophantine #approximations offered by NumerAppx to #sum(Cf(i)/pol(i),i=1..infinity) DeltaSeq:=proc(pol,Cf,n,L) local mu,lu,i: Digits:=100: lu:=[]: mu:=evalf(sum(subs(n=i,Cf/pol),i=1..infinity)): for i from degree(pol,n) to L do lu:=[op(lu),evalf(delt(NumerAppx(pol,Cf,n,i,i),mu),6)]: od: lu: end: #DeltaSeq1(pol,Cf,n,L,Const): The sequence of deltas for the diophantine #approximations offered by NumerAppx to the constant Const that #should be: sum(Cf(i)/pol(i),i=1..infinity) DeltaSeq1:=proc(pol,Cf,n,L,Const) local lu,i: Digits:=100: lu:=[]: for i from degree(pol,n) to L do lu:=[op(lu),evalf(delt(NumerAppx(pol,Cf,n,i,i),Const),6)]: od: lu: end: