#You can always re-set ERIC:=BIGGER;

ERIC:=300:

######################################################################
##AperyRecurrence: Save this file as AperyRecurrence.                #
#Must also have the package EKHAD in the same directory              #
#To use it,                                                          #
#stay in the                                                         #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read AperyRecurrence : <Enter>                     # 
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg@math.Rutgers.edu.                                          # 
#######################################################################


read  EKHAD:
 
 
print(` Version of May 14, 2002`):
print(`First version, ca. Sept. 4, 2001`):

print(`This is AperyRecurrence, 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.rutgers.edu/~zeilberg.`):
print(`Please report all bugs to: zeilberg@math.rutgers.edu .`):
lprint(``):
print(`Written by Doron Zeilberger`):
print(``):
print(`For general help type: Ezra();`):

Ezra1:=proc()
if args=NULL then
print(`This is AperyRecurrence, 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(``):
lprint(``):
print(` All the  Procedures are: `):
print(` ApT,beta, delta, DeltaOp, GenSeq, `):
print(` GenSeqEr,GenSeqP `):
print(` LcmM, Rigdelta, RogerGeneral, RogerS, RogerErS , Roger, RogerZ `):
print(` Roger1, SeqRat, SeqRatCC, TovAp`):
print(``):
print(` For specific help type "Ezra(procedure_name)" `):
else
Ezra(args):
fi:

end:
 
Ezra:=proc()
if args=NULL then
print(`This is AperyRecurrence, 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(``):
lprint(``):
print(` For specific help type "Ezra(procedure_name)"; For all procedures `):
print(` type: Ezra1(); `):
print(``):
print(` The MAIN Procedures are:`):  
print(`   Roger , RogerZ`):
print(``):
fi:


if nops([args])=1 and op(1,[args])=ApT then
print(`ApT(d,n,k): given a discrete function d(n,k) `):
print(`ApT(d,n,k) is the Apery-transform `):
fi:

if nops([args])=1 and op(1,[args])=beta then
print(`beta(ope1,N,n): Given a recurrence operator, ope1, in the`):
print(`shift operator N, and the variable n, finds the beta such`):
print(`that the sequences satisfying the recurrence are O((beta)^n)`):
fi:


if nops([args])=1 and op(1,[args])=delta then
print(`delta(sid): Given a sequence of rational numbers, sid`):
print(`finds the smallest delta's s.t. |K-an/bn|<=C/bn^(1+delta)`):
print(`where K is its estimated limit, i.e. its last term`):
print(`For example, try: with(combiant):`):
print(`delta([seq(fibonacci(i+1)/fibonacci(i),i=1..100)]):`):
fi:

if nops([args])=1 and op(1,[args])=DeltaOp then
print(`DeltaOp(ope,N):If a(n) and b(n) are two solutions`):
print(`of the second  order operator ope (in the discrete variable n`):
print(`and the shift N), finds the first order recurrence satisfied by`):
print(`Delta(n)=a(n+1)*b(n)-a(n)*b(n+1) `):
fi:

if nops([args])=1 and op(1,[args])=`GenSeq` then
print(`GenSeq(rec,n,N,ini,K): Given a recurrence operator, rec(n,N)`):
print(`in terms of the variable n and the shift operator N,`):
print(`and a list of initial values, ini, generates the sequence`):
print(`from n=0 to n=K`):
print(`For example, try: GenSeq(N^2-N-1,n,N,[0,1],10)`):; 
fi:


if nops([args])=1 and op(1,[args])=`GenSeqEr` then
print(`GenSeqEr(rec,n,N,ini1,ini2,Er,L): Given a recurrence operator, `):
print(`rec(n,N) in terms of the variable n and the shift operator N,`):
print(`and two lists of initial values, ini1, inii2, generates the `):
print(` sequences (let's call them a[i], b[i] from n=0 up to when `):
print(` a[n+1]/b[n+1]-a[n]/b[n] is less then Er*10^(-5)`):
print(` but with at least L terms`):
print(` For example, try: GenSeqEr(N^2-N-1,n,N,[1,0],[0,1],1/10^10,24); `):
fi:

if nops([args])=1 and op(1,[args])=GenSeqP then
print(`GenSeqP(rec,n,N,ini,K,P): Given a recurrence operator, rec(n,N)`):
print(`in terms of the variable n and the shift operator N,`):
print(`and a list of initial values, ini, generates the sequence`):
print(`times lcm(P(1),..., P(n)), where P is a polynomial in n`):
print(`from n=0 to n=K`):
print(`For example, try: GenSeqP(N^2-N-1,n,N,[0,1],10,n); `):
fi:

if nops([args])=1 and op(1,[args])=LcmM then
print(`LcmM(p,n,n0): Given a polynomial p, in n, finds`):
print(`lcm(p(1),p(2), ...p(n0))`):
fi:

if nops([args])=1 and op(1,[args])=`Rigdelta` then
print(`Rigdelta(F,k,n,c): Given a binomial coefficient sum (w.r.t. k)`):
print(`with summand F (that depends on n and k), finds the rigorous delta`):
print(`(irrationality measure) of the associated Apery constant`):
print(`assuming that the Lcm has growth exp(c*n)`):
print(`it returns 0 if the operator is not second-order`):
print(`For example try:Rigdelta(binomial(n+k,k)*binomial(n,k),k,n,1);`):
fi:

if nops([args])=1 and op(1,[args])=`RogerS` then
print(`RogerS(F,k,n,K) given a proper hypergeometric term`):
print(`F, in k,n, first checks whetehr (i) the sequence `):
print(`b(n):=sum(F(k,n),k=0..n)`):
prit(`satisfies a second order recurrence (ii) if its associated sequence`):
print(`with initial values [0,1] is called a(n), finds the `):
print(`estimated limit of the sequence of ratios a(n)/b(n), `):
print(` followed by the estimated minimal delta`):
print(`for the sequence a(n)/b(n) `):
print(`up to the term n=K`):
print(`it returns 0 if the recurrence is not 2nd order`):
print(`For example, try RogerS(binomial(n,k)*binomial(n+k,k),k,n,40);`):
fi:

if nops([args])=1 and op(1,[args])=`RogerErS` then
print(`RogerErS(F,k,n,low1,high1, Er,L) given a proper hypergeometric term`):
print(`F, in k,n, first checks whetehr (i) the sequence `):
print(` b(n):=sum(F(k,n),k=low1..high1)`):
print(` (where low1 and high1 are the natural limits)`):
print(`satisfies a second order recurrence (ii) if its associated sequence`):
print(`with initial values [0,1] is called a(n), finds the `):
print(`limit (up to the error Er, but with at least L terms) and the `):
print(`estimated minimal delta`):
print(`for the sequence a(n)/b(n) `):
print(`up to the prescribed error`):
print(`It returns 0 if the recurrence is not 2nd order`):
print(`For example, `):
print(`try RogerErS(binomial(n,k)*binomial(n+k,k),k,n,0,n,1/10^12,30);`):
fi:



if nops([args])=1 and op(1,[args])=Roger then
print(`Roger(F,k,n,Di,L) given a proper hypergeometric term`):
print(`F, in k,n, (with one factor being binomial(n,k)`):
print(`finds the recurrence for b(n)=sum(b(n,k),k)`):
print(`  using  zeil then`):
print(`finds the associated sequence`):
print(`with initial values [1,0,...] let's call it a(n), then finds the `):
print(`limit (up to the error 1/10^Di, and at least L terms) and the `):
print(`estimated minimal delta`):
print(`for the sequence a(n)/b(n) `):
print(`up to the prescribed error`):

print(`For example, try; Roger(binomial(n,k)^5,k,n,12,30);`):
fi:

if nops([args])=1 and op(1,[args])=RogerGeneral then
print(`RogerGeneral(F,k,n,Di,L,INI) given a proper hypergeometric term`):
print(`F, in k,n, (with one factor being binomial(n,k)`):
print(`finds the recurrence for b(n)=sum(b(n,k),k)`):
print(`  using  zeil then`):
print(`finds the associated sequence`):
print(`with initial values INI let's call it a(n), then finds the `):
print(`limit (up to the error 1/10^Di, and at least L terms) and the `):
print(`estimated minimal delta`):
print(`for the sequence a(n)/b(n) `):
print(`up to the prescribed error`):

print(`For example, try; RogerGeneral(binomial(n,k)^5,k,n,12,30,[1,0,0]);`):
fi:

if nops([args])=1 and op(1,[args])=RogerZ then
print(`RogerZ(F,z,Er,L): Given a Laurent polynomial, F, in the variable z, `):
print(` and an error, Er, and an integer L, finds an appx., followed by `):
print(`the empirical delta, for the limit of a(n)/b(n) as n goes to`):
print(` infinity, where b(n):=Coeff. of z^0 of F^n and a(n) is the`):
print(`sequence satisfying the same recurrence as b(n) but with initial`):
print(`conditions, [1,0,0,..0] `):
print(`For example, try: RogerZ(1/z+3+2*z,z,1/10^20,40); `):
fi:

if nops([args])=1 and op(1,[args])=Roger1 then
print(`Roger1(F,k,n,low1,high1,Er,L) given a proper hypergeometric term`):
print(`F, in k,n,finds the recurrence for b(n)=sum(b(n,k),k=low1,high1)`):
print(`  using  zeil then`):
print(`finds the associated sequence`):
print(`with initial values [1,0,...] let's call it a(n), then finds the `):
print(`limit (up to the error Er, and at least L terms) and the `):
print(`estimated minimal delta`):
print(`for the sequence a(n)/b(n) `):
print(`up to the prescribed error`):

print(`For example, try; Roger1(binomial(n,k)^5,k,n,0,n,1/10^12,30);`):
fi:

if nops([args])=1 and op(1,[args])=SeqRat then
print(`SeqRat(rec,n,N,K): Given a recurrence operator rec in n and N`):
print(`and an integer K (>=25) computes the sequence of rations a[n]/a[n-1]`):
print(`for the sequence annihilating rec(n,N) and such that`):
print(`a(i)=0 for i<ORDER-1 and a(ORDER-1)=1, followed by empirical delta, `):
print(`For example, try: SeqRat(N^2-N-1,n,N,N,30); `):
fi:

if nops([args])=1 and op(1,[args])=SeqRatCC then
print(`SeqRatCC(rec,N,K): Given a recurrence operator `):
print(`with constant coefficients, rec in n and N`):
print(`and an integer K (>=25) computes the sequence of rations a[n]/a[n-1]`):
print(`for the sequence annihilating rec(n,N) and such that`):
print(`a(i)=0 for i<ORDER-1 and a(ORDER-1)=1`):
print(`followed by empirical delta, followed by rigorous delta`):
print(`For example, try: SeqRatCC(N^2-N-1,N,40);`):
fi:

if nops([args])=1 and op(1,[args])=TovAp then
print(`TovAp(rec,n,N,ini,K,p): The sequence GenSeq(rec,n,N,ini,K)`):
print(`times Lcm(p(1),..., p(n))`):
fi:

end:

#delta(sid): Given a sequence of rational numbers, sid
#finds the smallest delta in the list s.t. |K-an/bn|<=C/bn^(1+delta)
#where K is its estiamted limit, the last term
#For example, try: with(combiant):
#delta([seq(fibonacci(i+1)/fibonacci(i),i=1..100)]):
delta:=proc(sid)
local K,n,mu,i:
Digits:=max(ERIC,Digits):

n:=nops(sid):
if n<22 then
ERROR(`the sequence needs at least 22 terms`):
fi:
K:=evalf(sid[n]):

mu:=[]:
for i from max(10,n-20) to n-10 do
mu:=[op(mu),evalf(-log(abs(K-sid[i]))/log(denom(sid[i])))-1]:
od:
min(op(mu)):
end:


#GenSeq(rec,n,N,ini,K): Given a recurrence operator, rec(n,N)
#in terms of the variable n and the shift operator N,
#and a list of initial values, ini, generates the sequence
#from n=0 to n=K
#For example, try: GenSeq(N^2-N-1,n,N,[0,1],10); 
GenSeq:=proc(rec,n,N,ini,K)
local ORDER,gu,ope,i,lu,n0:
ORDER:=degree(rec,N):
if ORDER<>nops(ini) then
ERROR(`The fourth argument should have`.ORDER. `terms`):
fi:
ope:=expand(subs(n=n-ORDER,rec)/N^ORDER):
gu:=ini:
for n0 from nops(ini) to K do
lu:=0:
for i from ldegree(ope,N) to -1 do
lu:=lu+subs(n=n0,coeff(ope,N,i))*gu[nops(gu)+i+1]:
od:
lu:=-lu/subs(n=n0,coeff(ope,N,0)):
gu:=[op(gu),lu]:
od:
gu:
end:

ka:=proc(mu,L)
local i,gu,cu:
gu:=trunc(mu):
cu:=mu-trunc(mu):
for i from 1 to L do
gu:=gu+trunc(cu*10)/10^i:
cu:=mu-trunc(cu*10)/10^i:
od:  
gu:
end:

#beta(ope1,N,n): Given a recurrence operator, ope1, in the
#shift operator N, and the variable n, finds the beta such
#that the sequences satisfying the recurrence are O((beta)^n)
beta:=proc(ope1,N,n)
local gu,deg,gu1,lu:
gu:=expand(ope1):
deg:=degree(gu,n):
gu1:=coeff(gu,n,deg):
lu:=solve(gu1,N):
lu:=evalf(lu):
lu:=max(lu):
end:


#LcmM(p,n,n0): Given a polynomial p, in n, finds
#lcm(p(1),p(2), ...p(n0))
LcmM:=proc(p,n,n0) local i:
lcm(seq(subs(n=i,p),i=1..n0)):
end:


#GenSeqP(rec,n,N,ini,K,P): Given a recurrence operator, rec(n,N)
#in terms of the variable n and the shift operator N,
#and a list of initial values, ini, generates the sequence
#times lcm(P(1),..., P(n)), where P is a polynomial in n
#from n=0 to n=K
#For example, try: GenSeqP(N^2-N-1,n,N,[0,1],10,n); 
GenSeqP:=proc(rec,n,N,ini,K,P)
local gu,mu,i:
gu:=GenSeq(rec,n,N,ini,K):
mu:=[seq(LcmM(P,n,i),i=0..K)]:
[seq(mu[i]*gu[i],i=1..K+1)]:
end:




#GenSeqEr(rec,n,N,ini1,ini2,Er,L): Given a recurrence operator, rec(n,N)
#in terms of the variable n and the shift operator N,
#and two lists of initial values, ini1, inii2, generates the sequences
#(let's call them a[i], b[] from n=0 up to when 
#a[n+1]/b[n+1]-a[n]/b[n] is less then Er*10^(-5) and at least L terms
#For example, try: GenSeqEr(N^2-N-1,n,N,[1,0],[0,1],1/10^10,24); 
GenSeqEr:=proc(rec,n,N,ini1,ini2,Er,L)
local ORDER,gu1,gu2,ope,i,tu,n0,lu1,lu2:
ORDER:=degree(rec,N):
if ORDER<>nops(ini1) or ORDER<>nops(ini2) then
ERROR(`The fourth argument should have`.ORDER. `terms `):
fi:
ope:=expand(subs(n=n-ORDER,rec)/N^ORDER):
gu1:=ini1:
gu2:=ini2:
tu:=1:

for n0 from nops(ini1) to L  do
lu1:=0:
lu2:=0:
for i from ldegree(ope,N) to -1 do
lu1:=lu1+subs(n=n0,coeff(ope,N,i))*gu1[nops(gu1)+i+1]:
lu2:=lu2+subs(n=n0,coeff(ope,N,i))*gu2[nops(gu2)+i+1]:
od:
lu1:=-lu1/subs(n=n0,coeff(ope,N,0)):
lu2:=-lu2/subs(n=n0,coeff(ope,N,0)):
gu1:=[op(gu1),lu1]:
gu2:=[op(gu2),lu2]:
tu:=gu1[nops(gu1)]/gu2[nops(gu2)]-gu1[nops(gu1)-1]/gu2[nops(gu2)-1]:
od:




for n0 from L+1 while abs(tu)>Er/10^5  do
lu1:=0:
lu2:=0:
for i from ldegree(ope,N) to -1 do
lu1:=lu1+subs(n=n0,coeff(ope,N,i))*gu1[nops(gu1)+i+1]:
lu2:=lu2+subs(n=n0,coeff(ope,N,i))*gu2[nops(gu2)+i+1]:
od:
lu1:=-lu1/subs(n=n0,coeff(ope,N,0)):
lu2:=-lu2/subs(n=n0,coeff(ope,N,0)):
gu1:=[op(gu1),lu1]:
gu2:=[op(gu2),lu2]:
tu:=gu1[nops(gu1)]/gu2[nops(gu2)]-gu1[nops(gu1)-1]/gu2[nops(gu2)-1]:
od:



gu1,gu2:
end:

#RogerErS(F,k,n,low1,high1,Er,L) given a proper hypergeometric term
#F, in k,n, first checks whetehr (i) the sequence b(n):=sum(F(k,n),k=0..n)
#satisfies a second order recurrence (ii) if its associated sequence
#with initial values [0,1] is called a(n), finds the 
#limit (up to the error Er, and at least L terms) and the 
#estimated minimal delta
#for the sequence a(n)/b(n) 
#up to the prescribed error
#it returns 0 if the recurrence is not 2nd order
#For example, try RogerErS(binomial(n,k)*binomial(n+k,k),k,n,1/10^12);
RogerErS:=proc(F,k,n,low1,high1,Er,L)
local lu,sid1,sid2,N,i,ini1,ini2,sid,Sid,K,mu,k1:
read `C:\\mekhkar\\irra\\EKHAD`:
lu:=zeil(F,k,n,N,2):
if lu=0 then
 RETURN(0):
fi:

lu:=lu[1]:

ini1:=[0,1]:

ini2:=[]:

mu:=0:
for k1 from subs(n=0,low1) to subs(n=0,high1) do
mu:=mu+eval(subs({n=0,k=k1},F)):
od:
ini2:=[op(ini2),mu]:

mu:=0:
for k1 from subs(n=1,low1) to subs(n=1,high1) do
mu:=mu+eval(subs({n=1,k=k1},F)):
od:
ini2:=[op(ini2),mu]:

Sid:=GenSeqEr(lu,n,N,ini1,ini2,Er,L):

sid1:=Sid[1]:
sid2:=Sid[2]:
K:=nops(sid1):

sid:=[seq(sid1[i]/sid2[i],i=1..K)]:
evalf(sid[K],10),evalf(delta(sid,K),10):
end:


#DeltaOp(ope,N):If a(n) and b(n) are two solutions
#of the second  order operator ope (in the discrete variable n
#and the shift N), finds the first order recurrence satisfied by
#Delta(n):=a(n+1)*b(n)-a(n)*b(n+1)
DeltaOp:=proc(ope,N)
local p,r:

if not degree(ope,N)=2 then
 ERROR(`Operator not second-order`):
fi:
p:=coeff(expand(ope),N,2):
r:=coeff(expand(ope),N,0):
p*N-r:
end:


 

#Rigdelta(F,k,n,c): Given a binomial coefficient sum (w.r.t. k)
#with summand F (that depends on n and k), finds the rigorous delta
#(irrationality measure) of the associated Apery constant
#assuming that the Lcm has growth exp(c*n)
#it returns 0 if the operator is not second-order
#For example try:Rigdelta(binomial(n+k,k)*binomial(n,k),k,n,1);

Rigdelta:=proc(F,k,n,c)
local ope,b,b1,N:
ope:=expand(zeil(F,k,n,N)[1]):
if degree(ope,N)<>2 then
 print(`Operator not second-order`):
 RETURN(0):
fi:

b:=evalf(beta(ope,N,n)):

b1:=evalf(abs(beta(DeltaOp(ope,N), N,n))):
evalf((log(b)-log(b1)-c)/(log(b)+c)):
end:


#TovAp(rec,n,N,ini,K,p): The sequence GenSeq(rec,n,N,ini,K)
#times Lcm(p(1),..., p(n))
TovAp:=proc(rec,n,N,ini,K,p)
local gu,n0:
gu:=GenSeq(rec,n,N,ini,K):

[seq(LcmM(p,n,n0)*gu[n0],n0=1..K)]:
end:


Roger1:=proc(F,k,n,low1,high1,Er,L)
local lu,sid1,sid2,N,i,ini1,ini2,sid,Sid,K,mu,k1,ORDER,i1,n1:
lu:=zeil(F,k,n,N):

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

lu:=lu[1]:

ORDER:=degree(lu,N):

if ORDER<2 then
RETURN(0):
fi:

ini1:=[seq(0,i1=1..ORDER-1),1]:

ini2:=[]:
for n1 from 0 to ORDER-1 do
mu:=0:
for k1 from subs(n=n1,low1) to subs(n=n1,high1) do
mu:=mu+eval(subs({n=n1,k=k1},F)):
od:
ini2:=[op(ini2),mu]:
od:

Sid:=GenSeqEr(lu,n,N,ini1,ini2,Er,L):

sid1:=Sid[1]:
sid2:=Sid[2]:
K:=nops(sid1):

sid:=[seq(sid1[i]/sid2[i],i=1..K)]:
evalf(sid[K]),delta(sid,K):
end:



#ApT(d,n,k): given a discrete function d(n,k) 
#ApT(d,n,k) is the Apery-transform 
ApT:=proc(d,n,k)
local gu,k1:
gu:=0:
for k1 from 0 to k do
gu:=gu+binomial(n,k1)*binomial(k,k1)*d(n,k1):
od:
gu:
end:


#SeqRat(rec,n,N,K): Given a recurrence operator rec in n and N
#and an integer K (>=25) computes the sequence of rations a[n]/a[n-1]
#for the sequence annihilating rec(n,N) and such that
#a(i)=0 for i<ORDER-1 and a(ORDER-1)=1
#followed by empirical delta, 
SeqRat:=proc(rec,n,N,K)
local lu,gu,i1,i:
if K<25 then
ERROR(`K>=25`):
fi:
Digits:=max(Digits,ERIC):
gu:=GenSeq(rec,n,N,[seq(0,i1=0..degree(rec,N)-2),1],K);
lu:=[seq(gu[i]/gu[i-1],i=degree(rec,N)+2..nops(gu))]:
lu,evalf(delta(lu),10):

end:

#SeqRatCC(rec,N,K): Given a recurrence operator 
#with constant coefficients, rec in n and N
#and an integer K (>=25) computes the sequence of rations a[n]/a[n-1]
#for the sequence annihilating rec(n,N) and such that
#a(i)=0 for i<ORDER-1 and a(ORDER-1)=1
#followed by empirical delta, followed by rigorous delta
#For example, try: SeqRatCC(N^2-N-1,N,40);
SeqRatCC:=proc(rec,N,K)
local mu,lu,gu,i1,n,i,mu1:
if K<25 then
ERROR(`K>=25`):
fi:
Digits:=max(ERIC,Digits):
gu:=GenSeq(rec,n,N,[seq(0,i1=0..degree(rec,N)-2),1],K);
lu:=[seq(gu[i]/gu[i-1],i=degree(rec,N)+2..nops(gu))]:

mu:=evalf([solve(rec,N)]):
mu1:=[seq(abs(mu[i1]),i1=1..nops(mu))]:

if abs(mu[1]-mu1[1])>1/10^30 then
print(`Largest root is not real`):
fi:
lu,evalf(delta(lu),10),evalf((log(mu1[1])-log(mu1[2]))/log(mu1[1])-1,10):

end:


Roger:=proc(F,k,n,Di,L)
local lu,sid1,sid2,N,i,ini1,ini2,sid,Sid,K,mu,k1,ORDER,i1,n1,Er:
Er:=1/10^(Di+1):
lu:=zeil(F,k,n,N):

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

lu:=lu[1]:

ORDER:=degree(lu,N):

if ORDER<2 then
RETURN(0):
fi:

ini1:=[seq(0,i1=1..ORDER-1),1]:

ini2:=[]:
for n1 from 0 to ORDER-1 do
mu:=0:
for k1 from 0 to n1 do
mu:=mu+eval(subs({n=n1,k=k1},F)):
od:
ini2:=[op(ini2),mu]:
od:

Sid:=GenSeqEr(lu,n,N,ini1,ini2,Er,L):

sid1:=Sid[1]:
sid2:=Sid[2]:
K:=nops(sid1):

sid:=[seq(sid1[i]/sid2[i],i=1..K)]:
evalf(sid[K],Di),evalf(delta(sid,K),5):
end:



RogerZ:=proc(F,z,Er,L)
local n,lu,sid1,sid2,N,i,ini1,ini2,sid,Sid,K,ORDER,i1,n1:
lu:=AZd(F^n/z,z,n,N):
if lu=0 then
 RETURN(0):
fi:

lu:=lu[1]:

ORDER:=degree(lu,N):

if ORDER<2 then
RETURN(0):
fi:

ini1:=[seq(0,i1=1..ORDER-1),1]:

ini2:=[]:
for n1 from 0 to ORDER-1 do
ini2:=[op(ini2),coeff(expand(F^n1),z,0)]:
od:

Sid:=GenSeqEr(lu,n,N,ini1,ini2,Er,L):

sid1:=Sid[1]:
sid2:=Sid[2]:
K:=nops(sid1):

sid:=[seq(sid1[i]/sid2[i],i=1..K)]:
evalf(sid[K]),evalf(delta(sid,K),6):
end:



#RogerGeneral(F,k,n,Di,L,INI): like RogerGeneral(F,k,n,Di,L) but with
#arbitrary initial conditions INI,
RogerGeneral:=proc(F,k,n,Di,L,ini1)
local lu,sid1,sid2,N,i,ini2,sid,Sid,K,mu,k1,ORDER,i1,n1,Er:
Er:=1/10^(Di+1):
lu:=zeil(F,k,n,N):

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

lu:=lu[1]:

ORDER:=degree(lu,N):

if ORDER<2 then
RETURN(0):
fi:



ini2:=[]:
for n1 from 0 to ORDER-1 do
mu:=0:
for k1 from 0 to n1 do
mu:=mu+eval(subs({n=n1,k=k1},F)):
od:
ini2:=[op(ini2),mu]:
od:

Sid:=GenSeqEr(lu,n,N,ini1,ini2,Er,L):

sid1:=Sid[1]:
sid2:=Sid[2]:
K:=nops(sid1):

sid:=[seq(sid1[i]/sid2[i],i=1..K)]:
evalf(sid[K],Di),evalf(delta(sid,K),5):
end: