#######################################################################
## MarkovWZdiag: save this file as MarkovWZdiag. To use it, stay in   #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read MarkovWZdiag : <Enter>                        #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Mohamud Mohammed and Doron Zeilberger,                  #
## Rutgers University ,                                               #
#######################################################################

Digits:=200:

#For log(2)
H1:=(-1)^z*z!/(x+z+1)!:

#For Zeta(2)
H2:=(-1)^z*2*z!/(x+z+1)/(2*x+z+1)!;

H2Old:=(z!/(x+z+1)!)^2:

#For Zeta(3)
H3:=z!/(2*x+z+1)!/(x+z+1)^2;
#For 1-Catalan
C1:=(-1)^z*(z+1/2)!^2/(2*x+z+3/2)!^2/4;
print(`First Version: May 28, 2004: tested for Maple 8 and 9 `):
print(`Version of May 28, 2004:  `):
print():
print(`This is MarkovWZdiag, A Maple package`):
print(`accompanying the article: The WZ-Markov Method `):
print(` by M. Mohammed and D. Zeilberger`):
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name)" `):
print():
print(): 
 
print(` For help on procedures related to the diagonal contour `):
print(`type "ezraD();"`):


ezra1:=proc()

if args=NULL then
print(`All the procedures are  `):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures: `):
print(`BestDelt, BestDeltAll, BestDeltT. BestKernelAppx,`):
print(` CheckMWZ, CheckMWZI `):
print(` CheckZeil, Delt, Gosper1 GosperP, GosperP1, Markov, `):
print(` SeqRat, SFR, VecRec, WtAveD, Zeil `):
else
ezra(args):
fi:
end:
 
 
ezraD:=proc()
if args=NULL then
print(`Contains procedures: `):
print(`MarkovD, AccDiag, AppxDiagf, AppxDiagr,`):
print( ):
print(`For a specific procedure type: ezraD(proc. name)`):

elif nargs=1 and args[1]=AppxDiagf then
print(`AppxDiagf(F,k,n,X0,N0): Given a closed-form F(n,k), outputs`):
print(`the summand of an infinite sum followed by a floating-point`):
print(`approximation using X0 walks, and N0 terms over`): 
print(`the diagonal contour.`):
print(` That means: sum(F(0,z),z=0..infinity)= sum(G(x,0),x=0..N0-1)+`):
print(`sum(G(N0+x,x)+F(N0+x+1,x),x=0..X0)`):
print(`For example, try: AppxDiagf((z!/(2*x+z+1)!)^2,z,x,10,10);`):

elif nargs=1 and args[1]=AppxDiagr then
print(`AppxDiagr(F,k,n,X0,N0): Given a closed-form F(n,k), outputs`):
print(`the summand of an infinite sum followed by a rational`):
print(`approximation using N0 walks, and X0 terms over`): 
print(`the diagonal contour.`):
print(` That means: sum(F(0,z),z=0..infinity)= sum(G(x,0),x=0..N0-1)+`):
print(`sum(G(N0+x,x)+F(N0+x+1,x),x=0..X0)`):
print(`For example, try: AppxDiagr((z!/(2*x+z+1)!)^2,z,x,10,10);`):


elif nargs=1 and args[1]=AccDiag then
print(` AccDiag(H,k,n,N0): Given a hypergemetric term H in discrete `):
print(`and an integer N0, the output consists of three parts`):
print(`The first part is the summand of the slowly converging-series`):
print(`that has to be accelerated, which is simply subs(x=0,H)`):
print(`the second part is the transition matrix (let's call it A)`):
print(`The third part is a list of length three,`):
print(`whose first component is coeff(RF(x,0),a_i),`):
print(`and the second component is  F(N0+x+1,x)`):
print(`the third component is coeff(RF(N0+x,x),a_i) ,`):
print(`For example, try: AccDiag((z!/(2*x+z+1)!)^2,z,x,10);`):
print(``):


elif nargs=1 and args[1]=MarkovD then
print(` MarkovD(H,k,n): Given a hypergemetric term H in discrete `):
print(` in x and z variables k and n, returns the WZ-Markov pair.`):
print(`The first component is a transition matrix A `):
print(`The second component is the coeff(RG,a_i)`):
print(`For example, try: MarkovD((z!/(2*x+z+1)!)^2,z,x);`):
 
else
 ezraD(args):
fi:
end:


ezra:=proc()

if args=NULL then
print(`MarkovWZ `):
print(`A Maple package for discovering Markov-WZ Pairs. `):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`For a list of all procedures, type ezra1(); `):
print(`The Main  procedures are: `):
print(`  ACC,`):
print(`Appx, Appxf, BestDeltT, cnkMWZ, delt, EAcc, EvalMWZ, Markov,`):
print(` WtAve, WtAve1, WtAve2, Zeil `):
print():

elif nargs=1 and args[1]=ACC then
print(` ACC(H,z,x): Given a hypergemetric term H in the discrete variables `):
print(` z and x, finds an acceleration scheme for sum(H(z,0),z=0..infinity`):
print(` implied by Markov(H,z,x) `):
print(`In other words, the output consists of three parts`):
print(`The first part is the summand of the slowly converging-series`):
print(`that has to be accelerated, which is simply subs(x=0,H)`):
print(`the second part is the transition matrix (let's call it A)`):
print(`defining the functions`):
print(`a_0(x), ..., a_{L-1}(x), such that`):
print(`[a_0(x+1), ... , a_{L-1}(x+1}]= A times `):
print(`[a_0(x), ... , a_{L-1}(x}] `):
print(`with the initial conditions [a_0(0),a_1(0), ..., a_{L-1}(0)(0)=`):
print(`[1,0,0,...,0] `):
print(`The third part is a vector of length L, lets' call it`):
print(`[b_0(x), ... , b_{L-1}(x)] such that the summand of the right`):
print(` hand side (the accelerating series) is `):
print(`a_0(x)b_0(x)+ ... +a_{L-1}(x) b_{L-1}(x) `):
print(`Warning: the Accelerating scheme is good only if the `):
print(`the hyptoheses of the theorem are met (that certain things `):
print(` go to 0, these hypotheses should be checked by the user`):
print(`For example, try ACC((z!/(x+z+1)!)^2,z,x); `):


elif nargs=1 and args[1]=Appx then
print(`Appx(F,k,n,N0): Given a closed-form F(n,k), outputs`):
print(`the summand of an infinite sum followed by a floating-point`):
print(`approximation using N0 terms. For example, try:`):
print(`Appx((z!/(x+z+1)!)^2,z,x,10);`):

elif nargs=1 and args[1]=Appxf then
print(`Appxf(H,k,n,L): Given a kernel H`):
print(`retuns subs(n=0,H) and `):
print(`the approximation, in floating, using the L terms of the`):
print(`fast series. For example, try:`):
print(`Appxf((z!/(x+z+1)!)^2,z,x);`):

elif nargs=1 and args[1]=BestDelt then
print(`BestDelt(F,C,k,n,r0): Given a discrete function F(k,n) such that `):
print(`sum(F(k,0),k=0..infinity)=C, finds the biggest `):
print(`of delt(cnkMWZ(MWZ(F,k,n),k,n,k0,n0),C): for k0+n0=r0 followed`):
print(`by the champion [k0,n0]`):
print(`For example, try: BestDelt(H1,log(2),z,x,20); `):
print(`or: BestDelt(H2,Zeta(2),z,x,20); `):
print(`or: BestDelt(H3,Zeta(3),z,x,20); `):
print(`or: BestDelt(C1,1-Catalan,z,x,20); `):

elif nargs=1 and args[1]=BestDeltAll then
print(`BestDeltAll(F,C,k,n,R0,R1): `):
print(`Given a discrete function F(k,n) such that `):
print(`sum(F(k,0),k=0..infinity)=C, finds `):
print(`BestDelt(F,C,k,n,r0) for r0=R0..R1`):
print(`For example, try: BestDeltAll(H1,log(2),z,x,10,20); `):
print(`or: BestDeltAll(H2,Zeta(2),z,x,10,20); `):
print(`or: BestDeltAll(H3,Zeta(3),z,x,10,20); `):
print(`or: BestDeltAll(C1,1-Catalan,z,x,10,20); `):

elif nargs=1 and args[1]=BestDeltT then
print(`BestDeltT(F,C,z,x,Pars,ParsLim,R0): Given a discrete function F(z,x)`):
print(`such that sum(F(z,0),z=0..infinity) is C, and F also`):
print(`depends on integer parameters in Pars, finds the`):
print(`best choice in the Pars betwen the limits given by ParsLim`):
print(`For example, try:`):
print(` BestDeltT((-1)^z*z!/(a*x+z+1)!,log(2),z,x,[a],[[1,3]],20)`):

elif nargs=1 and args[1]=BestKernelAppx then
print(`BestKernelAppx(r,z,x,N0): the kernel of the form`):
print(`H:=(z!/(2*x+z+1)!)^i/(x+z+1)^(r-i) (1<=i<=r) for`):
print(` approximating Zeta(r) that`):
print(`gives the best approximation after N0 terms, followed by the error`):
print(`For example, try: BestKernelAppx(4,z,x,40);`):

elif nargs=1 and args[1]=CheckMWZ then
print(`CheckMWZ(MWZ,k,n,L): Given a Markov-WZ pair MWZ in the`):
print(` variables k and n, and an integer L,  `):
print(`checks, numerically, for  range determined by L`):
print(`variables k,n, checks that sum(F(n,k),k=0..n)=1 for n<=L`):
print(`For example, try: CheckMWZ(Markov(binomial(n,k)^3,k,n),k,n,10)`):

elif nargs=1 and args[1]=CheckMWZI then
print(`CheckMWZI(MWZ,k,n,L): Given a Markov-WZ pair MWZ in the`):
print(`variables k,n, checks that sum(F(n,k),k=0..n)=1 for n<=L`):
print(`For example, try: CheckMWZI(Markov(binomial(n,k)^3,k,n),k,n,10)`):

elif nargs=1 and args[1]=CheckZeil then
print(` CheckZeil(F,k,n,N,ope,cert): checks the output to Zeil `);


elif nargs=1 and args[1]=cnkMWZ then
print(`cnkMWZ(MWZ,k,n,k0,n0): the Potential function c(n,k) for`):
print(`the Markov-WZ pair MWZ at k=k0 , n=n0`):

elif nargs=1 and args[1]=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)):`):

elif nargs=1 and args[1]=Delt then
print(`Delt(F,C,z,x,r0,s0): Given a discrete function F(z,x)`):
print(`in the variables z,x such that`):
print(`such that sum(F(z,0),z=0..infinity) is C, `):
print(`and F is a hypergeometric term, find the delta of`):
print(`cnk at r0,s0`):
print(`For example, try: Delt(H3,Zeta(3),z,x,10,10);`):

elif nargs=1 and args[1]=EAcc then
print(`EAcc(H,z,x): Explicit Acceleration Formula, if available`):
print(`otherwise returns FAIL. For example, try `):
print(`EAcc((-1)^z*z!/(x+z+1)!,z,x);`):

elif nargs=1 and args[1]=EvalMWZ then
print(`EvalMWZ(MWZ,k,n,k0,n0): Given a Markov-WZ pair`):
print(` MWZ=[H,Mat,Left,Right]`):
print(`in the variables k and n, evaluates its value (a list of lenght) 2`):
print(`at k=k0 and n=n0`):
print(`evaluates it at n=n0 and k=k0`):
print(`For example, try: EvalMWZ(Markov((z!/(x+z+1)!)^2,z,x),z,x,1,1);`):

elif nargs=1 and args[1]=Gosper1 then
print(`Gosper1(POL,F,k): Given a polynomial POL(k), and a hypergeometric`):
print(`F(k), and a variable k, decides whther there exist a hypergeometric`):
print(`z(k) such that z(k+1)-z(k)=POL(k)*F(k), or FAIL.`):
print(`For example, try: Gosper(1,k*k!,k)`):


elif nargs=1 and args[1]=GosperP then
print(` GosperP(POL,F,k,VAR1): Given a polynomial POL(k), `):
print(` featuring coefficients from a list VAR1, and a hypergeometric`):
print(` F(k), and a variable k, finds the general assignment for VAR1 `):
print(` that will be it Gosperable together with the the anti-difference `):
print(` z(k) divided by F `):
print(`z(k) such that z(k+1)-z(k)=POL(k)*F(k)`):
print(`For example, try: GosperP(a0+a1*k,k!,k,[a0,a1]); `):

elif nargs=1 and args[1]=GosperP1 then
print(` GosperP1(POL,F,k,VAR1): Like GosperP `):
print(`but only returns the solution set and the the z(k) `):
print(`For example, try: GosperP1(a0+a1*k,k!,k,[a0,a1]); `):


elif nargs=1 and args[1]=Markov then
print(` Markov(H,k,n): Given a hypergemetric term H in discrete `):
print(` variables k and n, returns the WZ-Markov pair.`):
print(`The first component if the kernel (the input H)`):
print(`The second component is a transition matrix A such that `):
print(` [a[0](n+1), a[1](n+1), ... ]= A times` ):
print(` [a[0](n),a[1](n), ... ] `):
print(`and [a[0](0), a[1](0), ...]=[1,0,0,...] `):
print(`The third component is the F part in the WZ-pair`):
print(`Given in terms of a vector to be dot-producted with`):
print(`[a[0](n),a[1](n),..], and the result multiplied by H`):
print(`The fourth component is the G part in the WZ-pair`):
print(`Given in terms of a vector to be dot-producted with`):
print(`[a[0](n),a[1](n),..], and the result multiplied by H`):
print(`i.e. we have F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k) `):
print(` For example, try: Markov(binomial(n,k)^3,k,n) `):

elif nargs=1 and args[1]=SeqRat then
print(`SeqRat(F,k,n,N0): Given a closed-form F, in variables k,n, `):
print(`gives the sequence of ratios in terms of Appx(F,k,n,N0) `):


elif nargs=1 and args[1]=SFR then
print(`SFR(A,n,ini,K): Given a list of lists A with functions of n`):
print(`a variable n and an initial vector such that vec(0)=ini`):
print(`finds vec(n) such that vec(n)=A(n)vec(n-1)`):
print(`Try SFR([[1]],n,10,[1]);`):

elif nargs=1 and args[1]=VecRec then
print(`VecRec(Mat,x,x0): Given a transition matrix in variable x`):
print(`computes the x0's term of the sequence A(x0+1)=M(x0)A(x0)`):
print(`where A(0)=[1,0,0, ..., 0];`):

elif nargs=1 and args[1]=WtAve then
print(`WtAve(MWZ,k,n,F,a0): the weighted average of the Potential funcion`):
print(`of the Markov-WZ pair MWZ (in variables k and n) and F `):
print(`over n+k=a0`):
print(`For example, try`):
print(`WtAve(Markov(H2,z,x),z,x,binomial(z+2*x,x)*binomial(x+z,z)^2,20);`):
print(`WtAve(Markov(H3,z,x),z,x,binomial(z+2*x,x)^2*binomial(x+z,z)^2,20);`):

elif nargs=1 and args[1]=WtAve1 then
print(`WtAve1(MWZ,k,n,F,a0): the weighted average of the Potential funcion`):
print(`of the Markov-WZ pair MWZ (in variables k and n) and F at a0`):
print(`over n=a0 0<=k<=a0`):
print(`For example, try`):
print(`WtAve1(Markov(H1,z,x),z,x,binomial(x,z)*binomial(x+z,z),20);`):

elif nargs=1 and args[1]=WtAve2 then
print(`WtAve2(MWZ,k,n,F,a0): the weighted average of the Potential funcion`):
print(`of the Markov-WZ pair MWZ (in variables k and n) and F at a0`):
print(`over k=a0 0<=dnk<=a0`):
print(`For example, try`):
print(`WtAve2(Markov(H1,z,x),z,x,binomial(x,z)*binomial(x+z,z),20);`):

elif nargs=1 and args[1]=WtAveD then
print(`WtAveD(H,k,n,F,a0): short for WtAve(Markov(H,k,n),k,n,F,a0);`):

elif nargs=1 and args[1]=Zeil then
print(` Zeil(F,k,n,N) : Given a proper hypergeometric term F in the`):
print(`discrete variables k and n, and a shift operator N `):
print(`returns the recurrence operator ope(N,n), and a G(n,k)/F(n,k)  `):
print(` such that ope(N,n)F(n,k)=G(n,k+1)-G(n,k) `):
print(`For example, try: Zeil(binomial(n,k),k,n,N); `):

else
 print(`There is no Ezra for: `, convert( args, string ) ):
fi:
 
end:







#Gosper1(POL,F,k): Given a polynomial POL(k), and a hypergeometric
#F(k), and a variable k, decides whther there exist a hypergeometric
#z(k) such that z(k+1)-z(k)=POL(k)*F(k)
Gosper1:=proc(POL,F,k)
local r,i,a,b,c,deg,x,eq,var,e,lu:
r:=normal(expand(simplify(subs(k=k+1,F)/F))):
a:=factor(numer(r)): b:=factor(denom(r)):
c:=POL:


lu:=MakeBest(a,b,c,k):
a:=lu[1]: b:=lu[2]:c:=lu[3]:


deg:=FindDeg(a,b,c,k):
if deg=FAIL then
 RETURN(FAIL):
fi:

x:=GenPol(k,deg,e):
var:=x[2]:
x:=x[1]:

lu:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x-c):
eq:={seq(coeff(lu,k,i),i=0..degree(lu,k))}:
var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

x:=subs(var,x):
normal(subs(k=k-1,b)*POL*x/c):
end:


#FindDeg(a,b,c,k): Finds the degree for x(k)
FindDeg:=proc(a,b,c,k)
local D1,D2,i1,A,B:
if degree(a,k)<>degree(b,k) or
   coeff(a,k,degree(a,k))<>coeff(b,k,degree(b,k)) then
    RETURN(degree(c,k)-max(degree(a,k),degree(b,k))):
fi:

A:=coeff(a,k,degree(a,k)-1):
B:=coeff(b,k,degree(b,k)-1):

D1:={degree(c,k)-degree(a,k)+1,(B-A)/coeff(a,k,degree(a,k))}:
D2:={}:
for i1 from 1 to nops(D1) do
 if type(D1[i1],integer) and D1[i1] >=0 then
    D2:=D2 union {D1[i1]}:
 fi:
od:
if D2={} then
  RETURN(FAIL):
else
 RETURN(max(op(D2))):
fi:
end:


#GenPol(n,deg,a): A generic polynomial of degree deg in n
#with coeffs. indexed by a, followed by the set of coeffs.
GenPol:=proc(n,deg,a) local i:
convert([seq(a[i]*n^i,i=0..deg)],`+`), {seq(a[i],i=0..deg)}: 
end:

#GosperP(POL,F,k,VAR1): Given a polynomial POL(k), 
#featuring coefficients from a list VAR1, and a hypergeometric
#F(k), and a variable k, finds the general assignment for VAR1
#that will be it Gosperable together with the the anti-difference
#z(k) divided by F
#z(k) such that z(k+1)-z(k)=POL(k)*F(k)
GosperP:=proc(POL,F,k,VAR1)
local r,i,a,b,c,deg,x,eq,var,e,lu,var1,tov,pu:
r:=normal(expand(simplify(subs(k=k+1,F)/F))):
a:=factor(numer(r)): b:=factor(denom(r)):

c:=POL:


lu:=MakeBest(a,b,c,k):
a:=lu[1]: b:=lu[2]:c:=lu[3]:

deg:=FindDeg(a,b,c,k):
if deg=FAIL then
 RETURN(FAIL):
fi:

x:=GenPol(k,deg,e):
var:=x[2]:
x:=x[1]:

lu:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x-c):
eq:={seq(coeff(lu,k,i),i=0..degree(lu,k))}:
var1:=solve(eq,var union convert(VAR1,set)):

if var1=NULL then
 RETURN(FAIL):
fi:

tov:={}:
for i from 1 to nops(var1) do
 if op(1,var1[i])=op(2,var1[i])    then
   tov:=tov union {op(1,var1[i])}:
  fi:
od:

x:=subs(var1,x):
pu:=[seq([VAR1[i],subs(var1,VAR1[i])],i=1..nops(VAR1))]:
normal(subs(k=k-1,b)*POL*x/c),pu, tov,var1:
end:


yafe:=proc(ope,N) local i:
sort(convert([seq(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N))],`+`)):end:

Zeil1:=proc(F,k,n,N,ORDER)
local gu,a,var,i,gu2,pu,POL,ope,cert:
var:=[seq(a[i],i=0..ORDER)]:
gu:=0:
for i from 0 to ORDER do
gu:=gu+a[i]*normal(expand(simplify(subs(n=n+i,F)/F))):
od:

POL:=numer(gu):
gu2:=denom(gu):

pu:=GosperP(POL,F/gu2,k,var):
if pu=FAIL then
 RETURN(0):
fi:
if pu[1]=0 then
 RETURN(0):
fi:

if nops(pu[3])>1 then
 print(`ORDER too big`):
fi:

ope:=convert([seq(factor(pu[2][i][2])*N^(i-1),i=1..ORDER+1)],`+`):
cert:=pu[1]/gu2:
ope:=subs(pu[3][1]=1,ope):
cert:=subs(pu[3][1]=1,cert):
ope:=normal(ope):
gu:=denom(ope):
cert:=normal(cert*gu):
ope:=numer(ope):
ope:=yafe(ope,N):
ope,cert:
end:

Zeil:=proc(F,k,n,N) local i,gu:
for i from 0 do
gu:=Zeil1(F,k,n,N,i):
if gu[1]<>0 then
 RETURN(gu):
fi:
od:
end:




Markov11:=proc(F,k,n,deg)
local gu,POL,POL1,F1,a,Var2, i,lu,i1,mu,mu1,cert,smol,yemin,pu,cert1:
option remember:

POL:=convert([seq(a[i](n)*k^i,i=0..deg)], `+`):
gu:=normal(normal(simplify(expand(subs(n=n+1,F)/F)))*subs(n=n+1,POL)-POL):
POL1:=numer(gu):
F1:=F/denom(gu):
Var2:=[seq(a[i](n+1),i=0..deg)]:

lu:=GosperP1(POL1,F1,k,Var2):
if lu=FAIL then
 RETURN(FAIL):
fi:
cert:=lu[2]/denom(gu):
lu:=lu[1]:
mu:=[]:

for i from 0 to deg do
mu1:=expand(subs(lu,a[i](n+1))):
mu:=[op(mu),[seq(factor(normal(coeff(mu1,a[i1](n),1))),i1=0..deg)]]:
od:
POL1:=subs(lu,POL):

smol:=normal(expand(subs(n=0,F))):
yemin:=subs(k=0,F):
cert1:=subs(k=0,cert):
pu:=[seq(coeff(cert1,a[i](n),1),i=0..deg)]:
mu,POL1,cert,[smol,yemin,pu]:
end:

#Gormim(P,k): all the roots a of P(k)=a such that in
#the fact
#it has non-linear factors
Gormim:=proc(P,k) local gu,mu,i,gu1:

if degree(P,k)=0 then
 RETURN({}):
fi:

mu:={}:
if degree(P,k)=1 then
RETURN({solve(P=0,k)}):
fi:

gu:=factor(P):

if type(gu,`*`) and degree(op(1,gu),k)=0 then
  gu:=factor(gu/op(1,gu)):
fi:

if type(gu,`^`) then
 gu1:=op(1,gu):
 if degree(gu1,k)<>1 then
   RETURN(FAIL):
 else
   RETURN({solve(gu1,k)}):
 fi:
fi:

for i from 1 to nops(gu) do
gu1:=op(i,gu):
 if type(gu1,`^`) then
  gu1:=op(1,gu1):
 fi:
 
 if degree(gu1,k)<>1 then
    RETURN(FAIL):
  else
    mu:=mu union {solve(gu1,k)}:
 fi:
od:
mu:
end:



Findh:=proc(a,b,k) local h,g,gu1,gu2,lu,i,j,lu1:

gu1:=Gormim(a,k):
gu2:=Gormim(b,k):
if gu1<>FAIL and gu2<>FAIL then
lu:={seq(seq(gu2[i]-gu1[j],j=1..nops(gu1)),i=1..nops(gu2))}:
 for i from 1 to nops(lu) do
  lu1:=op(i,lu):
   if type(lu1,integer) and lu1>0 then
     RETURN(lu1):
   fi:
 od:
  RETURN(FAIL):
else
print(`Using empirical approach`):
fi:
for h from 1 to 100 do
g:=gcd(a,subs(k=k+h,b)):
if degree(g,k)>0 then
RETURN(h):
fi:
od:
FAIL:
end:

MakeBetter:=proc(a,b,c,k) local h,g,a1,b1,c1,i1:
a1:=a: b1:=b:c1:=c:

h:=Findh(a,b,k):

if h<>FAIL then
g:=gcd(a1,subs(k=k+h,b1)):
 c1:=c1*convert([seq(subs(k=k-i1,g),i1=1..h)],`*`):
 a1:=normal(a1/g):
 b1:=normal(b1/subs(k=k-h,g)):
fi:

a1,b1,c1:
end:

MakeBest:=proc(a,b,c,k) local lu,lu1:
lu:=a,b,c:
lu1:=MakeBetter(lu,k):

while lu1<>lu do
lu:=lu1:
lu1:=MakeBetter(lu1,k):
od:
lu1:
end:





#GosperP1(POL,F,k,VAR1): Given a polynomial POL(k), 
#featuring coefficients from a list VAR1, and a hypergeometric
#F(k), and a variable k, finds the general assignment for VAR1
#that will be it Gosperable together with the the anti-difference
#z(k) divided by F
#z(k) such that z(k+1)-z(k)=POL(k)*F(k)
GosperP1:=proc(POL,F,k,VAR1)
local r,i,a,b,c,deg,x,eq,var,e,lu,var1:
r:=normal(expand(simplify(subs(k=k+1,F)/F))):
a:=factor(numer(r)): b:=factor(denom(r)):

c:=POL:


lu:=MakeBest(a,b,c,k):
a:=lu[1]: b:=lu[2]:c:=lu[3]:

deg:=FindDeg(a,b,c,k):
if deg=FAIL then
 RETURN(FAIL):
fi:

x:=GenPol(k,deg,e):
var:=x[2]:
x:=x[1]:

lu:=expand(a*subs(k=k+1,x)-subs(k=k-1,b)*x-c):
eq:={seq(coeff(lu,k,i),i=0..degree(lu,k))}:
var1:=solve(eq,var union convert(VAR1,set)):


if var1=NULL then
  RETURN(FAIL):
fi:

x:=subs(var1,x):


var1,factor(normal(subs(k=k-1,b)*POL*x/c)) :

end:




Markov1:=proc(F,k,n) local deg:
for deg from 0 while Markov11(F,k,n,deg)=FAIL do od:
Markov11(F,k,n,deg):end:


#CheckZeil(F,k,n,N,ope,cert): checks the output to Zeil

CheckZeil:=proc(F,k,n,N,ope,cert) local i:

normal(
convert([seq(coeff(ope,N,i)*simplify(
expand(subs(n=n+i,F)/F)),i=0..degree(ope,N))],
`+`)+
cert-subs(k=k+1,cert)*simplify(expand(subs(k=k+1,F)/F))):
end:


#SFR(A,n,ini,K): Given a list of lists A with functions of n
#a variable n and an initial vector such that vec(0)=ini
#finds vec(n) such that vec(n)=A(n)vec(n-1)
#Try SFR([[1]],n,10,[1]);
SFR:=proc(A,n,n0,ini)
local gu,A0,i,j:
option remember:
if n0=0 then
RETURN(ini):
fi:
gu:=SFR(A,n,n0-1,ini):
A0:=subs(n=n0-1,A):
[seq(convert([seq(A0[i][j]*gu[j],j=1..nops(A0[i]))],`+`),i=1..nops(ini))]:

end:



ACC:=proc(F,k,n) local gu,M1,M2,i,j:
option remember:
gu:=Markov1(F,k,n):
M1:=gu[1]:
M2:=gu[4][2]:
M2:=normal(expand(subs(n=n+1,M2)/M2)):

if abs(Appx1(F,k,n,20)[2]-Appx1(F,k,n,40)[2])>10^(-5) then
print(`The purorted acceleration scheme is a bad one`):
fi:

normal([gu[4][1],
[seq([seq(M1[i][j]*M2,j=1..nops(M1[i]))],i=1..nops(M1))],
simplify(subs({k=0,n=0},F))*gu[4][3]]):
end:

Acc:=proc(F,k,n) local gu,M1,M2,i,j:
option remember:

gu:=Markov1(F,k,n):
M1:=gu[1]:
M2:=gu[4][2]:
M2:=normal(expand(subs(n=n+1,M2)/M2)):

normal([gu[4][1],
[seq([seq(M1[i][j]*M2,j=1..nops(M1[i]))],i=1..nops(M1))],
simplify(subs({k=0,n=0},F))*gu[4][3]]):
end:


#VecRec(Mat,x,x0): Given a transition matrix in variable x
#computes the x0's term of the sequence A(x0+1)=M(x0)A(x0)
#where A(0)=[1,0,0, ..., 0];
VecRec:=proc(Mat,x,x0)
local k,i,j,lu,M1;
option remember:

k:=nops(Mat):
if x0=0 then
RETURN([1,0$(k-1)]):
fi:

lu:=VecRec(Mat,x,x0-1):
M1:=subs(x=x0-1,Mat):
[seq(convert([seq(M1[i][j]*lu[j],j=1..k)],`+`),i=1..k)]:

end:

VecRec1:=proc(Mat,x,x0,Vec1) local gu,lu,i:
gu:=VecRec(Mat,x,x0):
lu:=subs(x=x0,Vec1):
convert([seq(lu[i]*gu[i],i=1..nops(Vec1))],`+`):
end:

nthTerm:=proc(F,k,n,n0)
local lu:
lu:=Acc(F,k,n):
VecRec1(lu[2],n,n0,lu[3]):
end:

#Appx1(F,k,n,N0): Given a closed-form F(n,k), outputs
#the summand of an infinite sum followed by a floating-point
#approximation using N0 terms. For example, try:
#Appx((z!/(x+z+1)!)^2,z,x,10);
Appx1:=proc(F,k,n,N0) local n0:
normal(expand(subs(n=0,F))),
convert([seq(nthTerm(F,k,n,n0),n0=0..N0)],`+`):
end:

#Appx(F,k,n,N0): Given a closed-form F(n,k), outputs
#the summand of an infinite sum followed by a floating-point
#approximation using N0 terms. For example, try:
#Appx((z!/(x+z+1)!)^2,z,x,10);
Appx:=proc(F,k,n,N0) local n0,lu1,lu2:
lu1:=Appx1(F,k,n,N0):
lu2:=Appx1(F,k,n,trunc(N0/2)):
if abs(lu1-lu2)>10^(-5) then
ERROR(`Acc. scheme is no good, or N0 is too small`):
fi:

Appx1(F,k,n,N0):
end:

Appxf:=proc(F,k,n,N0) local lu:
lu:=Appx(F,k,n,N0): lu[1],evalf(lu[2]):end:

###Begin Elimination Part#####
LefMult:=proc(b,A) local i,j,k:
k:=nops(b):
if nops({seq(nops(A[i]),i=1..nops(A))})<>1 or nops(b)<>nops(A) then
ERROR(`Bad input)`):
fi:
[seq(expand(convert([seq(b[i]*A[i][j],i=1..k)],`+`)),j=1..nops(A[1]))]:
end:

Mult:=proc(B,A) local i:
[seq(LefMult(B[i],A) ,i=1..nops(B))]:
end:

#Eli(b,T,n,N): Given  a vector b of length k, say, of rational functions
#in  n, and a k by k matrix T with entries that are rational
#functions in n, finds the recurrence satisfied by
#x(n)=b(n).a(n) where a(n) is a k-vector solving
#a(n+1)=T(n)a(n)
Eli:=proc(b,T,n,N)
local i,k,A,a,gu,c,eq,var,ope:
k:=nops(b):
if nops(T)<>k or {seq(nops(T[i]),i=1..nops(T))}<>{k} then
ERROR(`Bad input`):
fi:

A:=[seq([a[i]],i=0..k-1)]:

gu:=[LefMult(b,A)]:


for i from 1 to k do
A:=Mult(subs(n=n+i-1,T),A):
gu:=[op(gu),LefMult(subs(n=n+i,b),A)]:
od:


var:={seq(c[i],i=0..k)}:
gu:=LefMult([seq(c[i],i=0..k)],gu)[1]:
eq:={seq(coeff(gu,a[i],1),i=0..k)}:
var:=solve(eq,var):
ope:=convert([seq(c[i]*N^i,i=0..k)],`+`):
ope:=subs(var,ope):

ope:=subs({seq(c[i]=1,i=0..k)},ope);
ope:=convert([seq(factor(coeff(ope,N,i))*N^i,i=0..k)],`+`):


end:

###End Elimination Part#####

Acc1:=proc(F,k,n,N) local lu:

lu:=Acc(F,k,n) :

lu[1],Eli(lu[3],lu[2],n,N):
end:

Bdok1:=proc(L,ope,n,N) local i,n0:
[seq(
convert([seq(subs(n=n0,coeff(ope,N,i))*L[n0+i+1],i=0..degree(ope,N))],`+`),
n0=0..nops(L)-degree(ope,N)-1)]:
end:


Bdok:=proc(F,k,n,N0)
local L,ope,N,n0:

L:=[seq(nthTerm(F,k,n,n0),n0=0..N0)]:
ope:=Acc1(F,k,n,N)[2]:
Bdok1(L,ope,n,N):
end:


##Floating point for compute

#AppxfOld(lu,k,n): Given an acceleration scheme lu
#returns the approximation with the prescribed Digits
#Appxf(Acc((z!/(x+z+1)!)^2,z,x),z,x);
AppxfOld:=proc(lu,k,n) local n0,su,te,eps:
eps:=10^(-Digits-3):

te:=nthTermf(lu,k,n,0):
su:=te:

for n0 from 1 while abs(te)>eps do
te:=nthTermf(lu,k,n,n0):
 su:=su+te:
od:
evalf(su,Digits-3):
end:


nthTermf:=proc(lu,k,n,n0):
VecRec1f(lu[2],n,n0,lu[3]):
end:


VecRec1f:=proc(Mat,x,x0,Vec1) local gu,lu,i:
gu:=VecRecf(Mat,x,x0):
lu:=subs(x=x0,Vec1):
convert([seq(lu[i]*gu[i],i=1..nops(Vec1))],`+`):
end:

#VecRecf(Mat,x,x0): Given a transition matrix in variable x
#computes the x0's term of the sequence A(x0+1)=M(x0)A(x0)
#where A(0)=[1,0,0, ..., 0];
VecRecf:=proc(Mat,x,x0)
local k,i,j,lu,M1;
option remember:

k:=nops(Mat):
if x0=0 then
RETURN([1.0,.0$(k-1)]):
fi:

lu:=VecRecf(Mat,x,x0-1):
M1:=subs(x=x0-1.0,Mat):
[seq(convert([seq(M1[i][j]*lu[j],j=1..k)],`+`),i=1..k)]:

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): 

if abs(evalf(c-a))<1/10^(Digits-3) or denom(a)<2 then
 RETURN(FAIL):
fi:
evalf(-log(abs(c-a))/log(denom(a)))-1: 
end:




MarkovA1:=proc(F,k,n,deg)
local gu,POL,POL1,F1,a,Var2, i,lu,i1,mu,mu1,cert:
option remember:

POL:=convert([seq(a[i](n)*k^i,i=0..deg)], `+`):
gu:=normal(normal(simplify(expand(subs(n=n+1,F)/F)))*subs(n=n+1,POL)-POL):
POL1:=numer(gu):
F1:=F/denom(gu):
Var2:=[seq(a[i](n+1),i=0..deg)]:

lu:=GosperP1(POL1,F1,k,Var2):
if lu=FAIL then
 RETURN(FAIL):
fi:
cert:=lu[2]/denom(gu):

lu:=lu[1]:
mu:=[]:

for i from 0 to deg do
mu1:=expand(subs(lu,a[i](n+1))):
mu:=[op(mu),[seq(factor(normal(coeff(mu1,a[i1](n),1))),i1=0..deg)]]:
od:
POL1:=subs(lu,POL):

POL1:=[seq(coeff(POL1,a[i1](n),1),i1=0..deg)]:
cert:=[seq(coeff(cert,a[i1](n),1),i1=0..deg)]:

[F,mu,POL1,cert]:

end:


Markov:=proc(F,k,n) local deg:
option remember:
for deg from 0 while MarkovA1(F,k,n,deg)=FAIL do od:
MarkovA1(F,k,n,deg):end:

#EvalMWZ(MWZ,k,n,k0,n0): Given a Markov-WZ pair MWZ=[H,Mat,Left,Right]
#in the variables k and n, evaluates its value (a list of lenght) 2
#at k=k0 and n=n0
#evaluates it at n=n0 and k=k0
EvalMWZ:=proc(MWZ,k,n,k0,n0) local i,H0,F0,G0,gu2,gu3,gu4:
gu2:=VecRec(MWZ[2],n,n0):
H0:=simplify(subs({n=n0,k=k0},MWZ[1])):
gu3:=simplify(subs({n=n0,k=k0},MWZ[3])):
gu4:=simplify(subs({n=n0,k=k0},MWZ[4])):
F0:=simplify(eval(H0*convert([seq(gu2[i]*gu3[i],i=1..nops(gu3))],`+`))):
G0:=simplify(eval(H0*convert([seq(gu2[i]*gu4[i],i=1..nops(gu3))],`+`))):
[F0,G0]:

end:

#CheckMWZ1(MWZ,k,n,k0,n0): Checks that the Markov-WZ pair MWZ
#is indeed a Markov-WZ pair at k=k0, n=n0
CheckMWZ1:=proc(MWZ,k,n,k0,n0)
evalb(EvalMWZ(MWZ,k,n,k0+1,n0)[2]-EvalMWZ(MWZ,k,n,k0,n0)[2]=
EvalMWZ(MWZ,k,n,k0,n0+1)[1]-EvalMWZ(MWZ,k,n,k0,n0)[1]):
end:


CheckMWZ:=proc(MWZ,k,n,L) local k0,n0:
{seq(seq(CheckMWZ1(MWZ,k,n,k0,n0),k0=L..L+20),n0=L+30..L+40)}:
end:


#cnkMWZ(MWZ,k,n,k0,n0): the Potential function c(n,k) for
#the Markov-WZ pair MWZ at k=k0 , n=n0
cnkMWZ:=proc(MWZ,k,n,k0,n0) local i:
add(EvalMWZ(MWZ,k,n,i,0)[1],i=0..k0-1)+
add(EvalMWZ(MWZ,k,n,k0,i)[2],i=0..n0-1):
end:


#cnkMWZ1(MWZ,k,n,k0,n0): the Potential function c(n,k) for
#the Markov-WZ pair MWZ at k=k0 , n=n0, both ways
cnkMWZ1:=proc(MWZ,k,n,k0,n0) local i:
add(EvalMWZ(MWZ,k,n,i,0)[1],i=0..k0)+
add(EvalMWZ(MWZ,k,n,k0+1,i)[2],i=0..n0),
add(EvalMWZ(MWZ,k,n,0,i)[2],i=0..n0)+
add(EvalMWZ(MWZ,k,n,i,n0+1)[1],i=0..k0):
end:


#WtAve(MWZ,k,n,F,a0): the weighted average of the Potential funcion5
#of the Markov-WZ pair MWZ (in variables k and n) and F at a0
#over k+n=a0
WtAve:=proc(MWZ,k,n,F,a0) local i:
eval(simplify(
convert(
[
seq(cnkMWZ(MWZ,k,n,i,a0-i)*subs({k=i,n=a0-i},F),i=0..a0)
],
`+`))
/
convert(
[
seq(simplify(subs({k=i,n=a0-i},F)),i=0..a0)
],
`+`)):
end:





#EAcc(H,z,x): Explicit Acceleration Formula, if available
#otherwise returns FAIL. For example, try EAcc((-1)^z*z!/(x+z+1)!,z,x);
EAcc:=proc(H,z,x)
local gu,M,f:
gu:=Markov(H,z,x):
M:=gu[2]:
if nops(M)<>1 or gu[3]<>[1] then
RETURN(FAIL):
fi:
M:=M[1][1]:
M:=rsolve({f(x+1)=M*f(x),f(0)=1},f(x)):
simplify([subs(x=0,gu[1]),subs(z=0,gu[1])*subs(z=0,gu[4][1])*M]):
end:

#WtAve1(MWZ,k,n,F,a0): the weighted average of the Potential funcion
#of the Markov-WZ pair MWZ (in variables k and n) and F at a0
#over n=a0 0<=k<=a0
#For example, try
#WtAve1(Markov(H1,z,x),z,x,binomial(x+1,z+1)*binomial(x+z+2,z+1),20);
WtAve1:=proc(MWZ,k,n,F,a0) local i:
simplify(
convert(
[
seq(cnkMWZ(MWZ,k,n,i,a0)*subs({k=i,n=a0},F),i=0..a0)
],
`+`)
/
convert(
[
seq(simplify(subs({k=i,n=a0},F)),i=0..a0)
],
`+`)):
end:

#WtAve2(MWZ,k,n,F,a0): the weighted average of the Potential funcion5
#of the Markov-WZ pair MWZ (in variables k and n) and F at a0
#over k=a0 0<=n<=a0
WtAve2:=proc(MWZ,k,n,F,a0) local i:
simplify(
convert(
[
seq(cnkMWZ(MWZ,k,n,i,a0)*subs({n=i,k=a0},F),i=0..a0)
],
`+`)
/
convert(
[
seq(simplify(subs({n=i,k=a0},F)),i=0..a0)
],
`+`)):
end:


#BestDelt(MWZ,C,k,n,r0): Given a Markov-WZ pair, MWZ, that is
#supposed to appx. a constant C, finds the biggest 
#of delt(cnkMWZ(MWZ,k,n,k0,n0),C): for k0+n0=r0 followed
#by the champion [k0,n0]
BestDelt:=proc(F,C,k,n,r0) local MWZ, aluf,shi, k0,muam:
MWZ:=Markov(F,k,n):
aluf:=[trunc((r0-1)/2),r0-trunc((r0-1)/2)]:
shi:=delt(cnkMWZ(MWZ,k,n,trunc((r0-1)/2),r0-trunc((r0-1)/2)),C):
if shi=FAIL then
 RETURN(FAIL):
fi:

for k0 from 1 to r0 do
 muam:=delt(cnkMWZ(MWZ,k,n,k0,r0-k0),C):
 if muma<>FAIL then
 if muam>shi then
   aluf:=[k0,r0-k0]:
   shi:=muam:
 fi:
 fi:
od:
shi,aluf:
end:



BestDeltAll:=proc(F,C,k,n,R0,R1) local r0:
[seq([BestDelt(F,C,k,n,r0)],r0=R0..R1)]:
end:


#CartProd(ParsLim): The Cartesian product of ParLim[1] x ParsLim[2]x ...
CartProd:=proc(ParsLim)
local gu,ParsLim1,k,akha,i,j:

k:=nops(ParsLim):
if k=1 then
 RETURN({seq([i],i=ParsLim[1][1]..ParsLim[1][2])}):
fi:
ParsLim1:=[op(1..k-1,ParsLim)]:
gu:=CartProd(ParsLim1):
akha:=ParsLim[k]:
{seq(seq([op(gu[i]),j],j=akha[1]..akha[2]),i=1..nops(gu))}:

end:

#BestDeltT(F,C,z,x,Pars,ParsLim,R0): Given a discrete function F(z,x)
#such that sum(F(z,0),z=0..infinity) is C, and F also
#depends on integer parameters in Pars, finds the
#best choice in the Pars betwen the limits given by ParsLim
#For example, try: BestDeltT((-1)^z*z!/(a*x+z+1)!,log(2),k,n,[a],[[1,3]],20)
BestDeltT:=proc(F,C,k,n,Pars,ParsLim,R0)
local gu,i,j,F1,lu,mu,shi,aluf,ma:
gu:=CartProd(ParsLim):
F1:=subs({seq(Pars[i]=gu[1][i],i=1..nops(gu[1]))},F):

lu:=BestDelt(F1,C,k,n,R0):
shi:=lu[1]:
aluf:=F1:
ma:=lu[2]:

for j from 2 to nops(gu) do
  F1:=subs({seq(Pars[i]=gu[j][i],i=1..nops(gu[j]))},F):
  mu:=BestDelt(F1,C,k,n,R0):
 if mu[1]>shi then
  aluf:=F1:
  ma:=mu[2]:
  shi:=mu[1]:
 fi:
od:
 
aluf,shi,ma:

end:


Delt:=proc(F,C,k,n,k0,s0):
delt(cnkMWZ(Markov(F,k,n),k,n,k0,s0),C):
end:


CheckMWZI:=proc(MWZ,k,n,L) local n0,k0:
{seq(add(EvalMWZ(lu,k,n,k0,n0)[1],k0=0..n0),n0=0..L)};
end:


#WtAveD(H,k,n,F,a0): short for WtAve(Markov(H,k,n),k,n,F,a0);
WtAveD:=proc(H,k,n,F,a0):
WtAve(Markov(H,k,n),k,n,F,a0):end:



#BestKernelAppx(r,z,x,N0): the kernel of the form
#H:=(z!/(2*x+z+1)!)^i/(x+z+1)^(r-i) (1<=i<=r) for approximating Zeta(r) that
#gives the best approximation after N0 terms, followed by the error
#For example, try: BestKernelAppx(4,z,x,40);

BestKernelAppx:=proc(r,z,x,N0) local H,i,shi,aluf,muam:
aluf:=1:
i:=1:
H:=(z!/(2*x+z+1)!)^i/(x+z+1)^(r-i);
shi:=evalf(Appxf(H,z,x,N0)[2]-Zeta(r)):

for i from 2 to r do
H:=(z!/(2*x+z+1)!)^i/(x+z+1)^(r-i);
muam:=evalf(Appxf(H,z,x,N0)[2]-Zeta(r)):
if muam<shi then
  aluf:=i:
  shi:=muam:
fi:
od:
(z!/(2*x+z+1)!)^aluf/(x+z+1)^(r-aluf),shi;
end:


#############  On Diagonal contour ################

MarkovD:=proc(F,k,n) local deg:
for deg from 0 while MarkovDiag(F,k,n,deg)=FAIL do od:
MarkovDiag(F,k,n,deg):end:

MarkovDiag:=proc(F,k,n,deg)
local gu,POL,POL1,F1,a,Var2, i,lu,i1,mu,mu1,cert,pu:
option remember:

POL:=convert([seq(a[i](n)*k^i,i=0..deg)], `+`):
gu:=normal(normal(simplify(expand(subs(n=n+1,F)/F)))*subs(n=n+1,POL)-POL):
POL1:=numer(gu):
F1:=F/denom(gu):
Var2:=[seq(a[i](n+1),i=0..deg)]:

lu:=GosperP1(POL1,F1,k,Var2):
if lu=FAIL then
 RETURN(FAIL):
fi:
cert:=lu[2]/denom(gu):
lu:=lu[1]:
mu:=[]:

for i from 0 to deg do
mu1:=expand(subs(lu,a[i](n+1))):
mu:=[op(mu),[seq(factor(normal(coeff(mu1,a[i1](n),1))),i1=0..deg)]]:
od:
POL1:=subs(lu,POL):


pu:=[seq(coeff(cert*F,a[i](n),1),i=0..deg)]:

mu,pu:

end:


# Returns G(x0,0)*vec((x0)
# which is the sum over the diagonal contour.

VecRecG:=proc(Mat,x,x0,Vec1) local gu,lu:

gu:=VecRec(Mat,x,x0):  

lu:=subs(x=x0,Vec1):

convert([seq(gu[i]*lu[i],i=1..nops(Vec1))],`+`):

end:


#Returns the n0th term of G(n0,0)

nthTermG:=proc(F,k,n,n0,N0)
local lu:
lu:=AccDiag(F,k,n,N0):

VecRecG(lu[2],n,n0,lu[3][1]):

end:


#AccDiag(F,k,n,N0):  returns F(0,z), A, 
#[coeff(FR(x,0),a_i), coeff(G(N0+x,x),a_i) , F(N0+x+1,x]

AccDiag:=proc(F,k,n,N0) local gu,M1,G0,mu,F1,G1,F0:
option remember:

gu:=MarkovD(F,k,n):

F0:=simplify(subs(n=0,F)):
mu:=gu[1]:
G1:=gu[2]:

G0:=subs(k=0,G1):

F1:=subs(n=N0+1+n,F):
F1:=subs(k=n,F1):

G1:=subs(n=N0+n,G1):
G1:=subs(k=n,G1):

normal(F0), mu, [normal(simplify(G0)),normal(F1),normal(G1)]: 

end:


# Returns F(N0+1+x0,x0)*vec((N0+1+x0) + RF(N0+x0,x0)*Vec(N0+x0)
# which is the sum over the diagonal contour.

VecRecDiag:=proc(Mat,x,x0,N0,F1,Vec1) local gu1,luG,luF,gu2,i:

gu1:=VecRec(Mat,x,N0+1+x0): # for the F(x0+1,x0)*vec((x0+1) part
luF:=subs(x=x0,F1):

gu2:=VecRec(Mat,x,N0+x0):   # for the RF(N0+x0,x0)*Vec(N0+x0) part
luG:=subs(x=x0,Vec1):

convert([seq(gu1[i]*luF,i=1..nops(Vec1))],`+`)+
convert([seq(gu2[i]*luG[i],i=1..nops(Vec1))],`+`):

end:


#Returns the n0th term of F(N0+n0+1,n0)+G(N0+n0,n0)

nthTermDiag:=proc(F,k,n,n0,N0)
local lu:

lu:=AccDiag(F,k,n,N0):

VecRecDiag(lu[2],n,n0,N0,lu[3][2],lu[3][3]):

end:


#AppxDiagr(F,k,n,N0,X0): Given a closed-form F(n,k), outputs
#the summand of an infinite sum followed by a rational
#approximation using \sum(G(x,0),x=0..N0-1) and X0 terms 
#over the diagonal contour.
#For example, try: AppxDiagr((z!/(x+z+1)!)^2,z,x,2,10);

AppxDiagr:=proc(F,k,n,N0,X0) local n0:
normal(expand(subs(n=0,F))),eval(
convert([seq(nthTermG(F,k,n,n0,N0),n0=0..N0-1)],`+`)),
eval(convert([seq(nthTermDiag(F,k,n,n0,N0),n0=0..X0)],`+`)):

end:


#AppxDiagf(F,k,n,N0,X0): Given a closed-form F(n,k), outputs
#the summand of an infinite sum followed by a floating-point
#approximation using \sum(G(x,0),x=0..N0-1) and X0 terms 
#over the diagonal contour.
#For example, try: AppxDiagf((z!/(x+z+1)!)^2,z,x,2,10);

AppxDiagf:=proc(F,k,n,N0,X0) local n0:
normal(expand(subs(n=0,F))),evalf(
convert([seq(nthTermG(F,k,n,n0,N0),n0=0..N0-1)],`+`)),
evalf(convert([seq(nthTermDiag(F,k,n,n0,N0),n0=0..X0)],`+`)):

end: