####################################################################### ## MultiAlmkvistZeilberger Save this file as MultiAlmkvistZeilberger # ## same directory, get into Maple (by typing: maple ) # ## and then type: read MultiAlmkvistZeilberger # ## Then follow the instructions given there # ## # ## Written by Mohamud Mohammed and Doron Zeilberger, # #Rutgers University , # ## [mohamudm, zeilberg] at math dot rutgers dot edu # ####################################################################### print(`First Written: Dec. 7, 2004: tested for Maple 8 and 9 `): print(`Version of Feb. 24, 2012 (adding verbose version) `): print(): print(`This is MultiAlmkvistZeilberger, one of the Maple packages`): print(`accompanying the article `): print(`Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and`): print(`the Sharpening of Wilf-Zeilberger Theory`): print(`by Mohamud Mohammed [not Moa Apagodu] and Doron Zeilberger `): print(``): 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(`For general help, and a list of the available functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): ezra:=proc() if nargs=0 then print(`MultiAlmkvistZeilberger`): print(`A Maple package for finding recurrences for multi-integrals of `): print(` Hypergeometric/Hyperexponential Integrands.`): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: DIAGpaper, DIAGpaperSR, MAZ, MAZpaper `): elif nargs=1 and args[1]=DIAGpaper then print(`DIAGpaper(TOP,BOT,a,x,n,pars,F,A,R): Inputs polynomials in the list of variables`): print(`x, TOP, and BOT, and a numeric or symbolic number a, and a symbol m, and a set of parameters pars`): print(`(if you are lazy, you can make it {}) and letters F,A,R, (for formulating the proof)`): print(`outputs a mathematical article`): print(`authored by Shalosh B. Ekhad that finds and proves a linear recurrence equation for the`): print(`diagonal sequence (called A(n)), i.e. the coefficient of x[1]^n*x[2]^n* in the formal power series`): print(`TOP/BOT^a. BOT must have a non-zero constant terms (for the formal power series to be well-defined)`): print(`For example, for Miklos Bona's problem, try:`): print(`DIAGpaper(1,(1-x)^2+(1-y)^2-1,1/2,[x,y],n,{},F,A,R):`): elif nargs=1 and args[1]=DIAGpaperSR then print(`DIAGpaperSR(TOP,BOT,a,x,n,A,n0,MaxC): Inputs polynomials in the list of variables`): print(`x, TOP, and BOT, and a numeric or symbolic number a, and a symbol m, and a set of parameters pars`): print(`and a letter A and pos. integers n0 and MaxC (for the guessing, if it returns FAIL makes them larger) `): print(`outputs a mathematical article`): print(`authored by Shalosh B. Ekhad that conjectures and semi-rigorously proves a linear recurrence equation for the`): print(`diagonal sequence (called A(n)), i.e. the coefficient of x[1]^n*x[2]^n* in the formal power series`): print(`TOP/BOT^a. BOT must have a non-zero constant terms (for the formal power series to be well-defined)`): print(`For example, for Miklos Bona's problem, try:`): print(`DIAGpaperSR(1,(1-x)^2+(1-y)^2-1,1/2,[x,y],n,A,30,10);`): elif nargs=1 and args[1]=MAZ then print(`MAZ(POL,H,rat,x,n,N,pars): Given a polynomial, POL, an hyper-exponential function, H, and`): print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`): print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`): print(`ouputs a recurrence operator, let's call it ope(n,N), and a multi-certificate, let's call it R`): print(`such that the function F(n;x):=POL*H*rat^n satisfies`): print(` ope(n,N)F=sum(D_{x[i]} (R[i]H*rat^n),i=1..nops(R)`): print(`In particular, try:`): print(`MAZ(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,A,{b,c});`): print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{});`): elif nargs=1 and args[1]=MAZpaper then print(`MAZpaper(POL,H,rat,x,n,pars,A,R): Given a polynomial, POL, an hyper-exponential function, H, and`): print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`): print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`): print(`and letters A and R for expressing the statement of the theorem`): print(`and proof`): print(`ouputs a paper stating and proving `): print(`A recurrence equation for the (contour or vanishing-at-endpoints) integral`): print(`of F:=POL*H*rat^n `): print(`with respect to the variables in the list [x] `): print(`Try:`): print(`MAZpaper(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,{b,c},A,R);`): print(`MAZpaper(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,{},A,R);`): else print(`There is no help for`, args): fi: end: ###Findrec #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then RETURN(FAIL): fi: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 0 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #end Findrec with(linalg): #findrecEx(f,DEGREE,ORDER,m1): Explores whether thre #is a good chance that there is a recurrence of degree DEGREE #and order ORDER, using the prime m1 #For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003); findrecEx:=proc(f,DEGREE,ORDER,m1) local ope,var,eq,i,j,n0,eq1,a,A1, D1,E1,Eq,Var,f1,n,N: option remember: f1:=f mod m1: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1: od: eq:= eq union {eq1}: od: Eq:= convert(eq,list): Var:= convert(var,list): D1:=nops(Var): E1:=nops(Eq): if E10 then RETURN(false): fi: if E1-D1>=1 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+1],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: if E1-D1>=2 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+2],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: true: end: ###End Findrec ez:=proc():print(`MAZ1(POL,H,rat,x,n,N,L,pars), MAZ(POL,H,rat,x,n,N,pars)`): print(`FindDeg(POL,H,rat,x,n,L)`): print(`CheckMAZ1(POL,H,rat,x,n,N,L,pars)`): end: _EnvAllSolutions:=true: ####From MultiZeilberger #IV1(d,n): all the vectors of non-negative integres of length d whose #sum is n IV1:=proc(d,n) local gu,i,gu1,i1: if d=0 then if n=0 then RETURN({[]}): else RETURN({}): fi: fi: gu:={}: for i from 0 to n do gu1:=IV1(d-1,n-i): gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}: od: gu: end: yafe:=proc(ope,N) local i: add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):end: #IV(d,n): all the vectors of non-negative integres of length d whose #sum is <=n IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end: #GenPol(kList,a,deg): The generic polynomial in kList of #degree deg, using the indexed variable a, followed by the set #of coeffs. GenPol:=proc(kList,a,deg) local gu,i,i1: gu:=IV(nops(kList),deg): convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`), i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}: end: #Extract1(POL,kList): extracts all the coeffs. of a POL in kList Extract1:=proc(POL,kList) local k1,kList1,i: k1:=kList[nops(kList)]: kList1:=[op(1..nops(kList)-1,kList)]: if nops(kList)=1 then RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}): else RETURN({seq( op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}): fi: end: ###########End from MultiZeilberger FindDeg:=proc(POL,H,rat,x,xSet,n,L) local s,t,h,e,i,Hbar,q,r,gu: s:=numer(rat): t:=denom(rat): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=normal(simplify(diff(Hbar,x)/Hbar)): q:=numer(gu): r:=denom(gu): max(degree(h,xSet)-degree( diff(r,x)+q,xSet)+degree(r,xSet) ,degree(h,xSet)-degree(r,xSet)+1+degree(r,xSet)): end: #MAZ1(POL,H,rat,x,n,N,L,pars): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1:=proc(POL,H,rat,x,n,N,L,pars) local deg,deg1,m,gu,i,i1: deg:=[]: for i from 1 to nops(x) do deg:=[op(deg),FindDeg(POL,H,rat,x[i],{seq(x[i1],i1=1..nops(x))},n,L)]: od: m:=min(op(deg)): for i from 0 to m do deg1:=[seq(deg[i1]-m+i,i1=1..nops(deg))]: gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,deg1): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: MAZ:=proc(POL,H,rat,x,n,N,pars) local gu,L: for L from 0 do gu:=MAZ1(POL,H,rat,x,n,N,L,pars): if gu<>FAIL then RETURN(gu): fi: od: end: CheckMAZ:=proc(POL,H,rat,x,n,N,pars) local F,gu,luL,luR,R,ope,i: F:=H*rat^n: gu:=MAZ(POL,H,rat,x,n,N,pars): if gu=FAIL then RETURN(FAIL): fi: ope:=gu[1]: R:=gu[2]: luL:=normal(add(subs(n=n+i,POL)*coeff(ope,N,i)*normal(simplify(subs(n=n+i,F)/F)),i=0..degree(ope,N))): luR:=normal(add( normal(diff(R[i]*F,x[i])/F ), i=1..nops(R))): evalb(normal(luL/luR)=1): end: #MAZ1deg(POL,H,rat,x,n,N,L,pars,deg): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1deg:=proc(POL,H,rat,x,n,N,L,pars,deg) local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,ku,i1,eqN,var1N,opeN,bilti,meka: s:=numer(rat): t:=denom(rat): ope:=add(e[i]*N^i,i=0..L): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=[]: for i from 1 to nops(x) do gu:=[op(gu),normal(simplify(diff(Hbar,x[i])/Hbar))]: od: q:=[]: r:=[]: for i from 1 to nops(gu) do q:=[op(q),numer(gu[i])]: r:=[op(r),denom(gu[i])]: od: var:={}: X:=[]: for i from 1 to nops(x) do ku:=GenPol(x,a[i],deg[i]): X:=[op(X),ku[1]]: var:=var union ku[2]: od: var:=var union {seq(e[i],i=0..L)}: gu:=h: for i from 1 to nops(x) do gu:=expand(normal(gu-(diff(r[i],x[i])+q[i])*X[i]/r[i]-r[i]*diff(X[i]/r[i],x[i]))): od: eq:=Extract1(numer(gu),x): eqN:=subs(n=9/17,eq): eqN:=subs({seq(pars[i]=1/(i^2+3),i=1..nops(pars))},eqN): var1N:=solve(eqN,var): opeN:=subs(var1N,ope): if opeN=0 then RETURN(FAIL): fi: var1:=solve(eq,var): ope:=subs(var1,ope): X:=subs(var1,X): if ope=0 then RETURN(FAIL): fi: bilti:={}: for i from 1 to nops(var1) do if op(1,op(i,var1))=op(2,op(i,var1)) then bilti:=bilti union {op(1,op(i,var1))}: fi: od: for i from 1 to nops(bilti) do if coeff(ope,bilti[i],1)<>0 then ope:=coeff(ope,bilti[i],1): X:=[seq(coeff(X[i1],bilti[i],1),i1=1..nops(X))]: meka:=coeff(ope,N,degree(ope,N)): ope:=expand(normal(ope/meka)): ope:=yafe(ope,N): X:=[seq(normal(X[i1]/meka/t^L),i1=1..nops(X))]: RETURN(ope,X): fi: od: end: #MAZ(POL,H,rat,x,n,pars,F,A,R) MAZpaper:=proc(POL,H,rat,x,n,pars,A,R) local gu,F,F1,ope,i,N: print(``): F1:=H*rat^n: gu:=MAZ(POL,H,rat,x,n,N,pars): ope:=gu[1]: if gu=FAIL then RETURN(FAIL): fi: print(`A Linear Recurrence Equation for an Integral in`, nops(x), ` variables `): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Let `): print(F(n,op(x))=F1): print(`and let`, A(n), `be its definite integral of`): print(POL*F(n,op(x))): print(`with respect to the continuous`): print(` variables `, op(x)): print(``): print(A(n) , `satisfies the following homogeneous linear recurrence with`): print(`polynomial coefficients: `): print(``): print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0): print(``): print(`Proof: We cleverly constuct the rational functions`): for i from 1 to nops(gu[2]) do print(R[i](n,op(x))=gu[2][i]): od: print(`With the motive that`): print(add(coeff(ope,N,i)*F(n+i,op(x)),i=0..degree(ope,N))= add(D[x[i]](R[i](n,op(x))*F(n,op(x))),i=1..nops(x))): print(``): print(`(Check!)`): print(`and the theorem follows upon integrating with respect to`, op(x)): print(` QED! `): end: #DIAGpaper(TOP,BOT,a,x,n,pars,F,A,R): Inputs polynomials in the list of variables #x, TOP, and BOT, and a numeric or symbolic number a, and a symbol m, and a set of parameters pars #(if you are lazy, you can make it {}) and letters F,A,R, (for formulating the proof) #outputs a mathematical article #authored by Shalosh B. Ekhad that finds and proves a linear recurrence equation for the #diagonal sequence (called A(n)), i.e. the coefficient of x[1]^n*x[2]^n* in the formal power series #TOP/BOT^a. BOT must have a non-zero constant terms (for the formal power series to be well-defined) #For example, for Miklos Bona's problem, try: #DIAGpaper(1,(1-x)^2+(1-y)^2-1,1/2,[x,y],n,{},F,A,R): DIAGpaper:=proc(TOP,BOT,a,x,n,pars,F,A,R) local gu,i,j,ope,N,t0: t0:=time(): if subs({seq(x[i]=0,i=1..nops(x))},BOT)=0 then print(`The bottom`, BOT, `has zero constant-terms , so the formal power series is undefined`): RETURN(FAIL): fi: gu:=MAZ(TOP,1/BOT^a/mul(x[i],i=1..nops(x)),1/mul(x[i],i=1..nops(x)),x,n,N,pars): print(``): ope:=gu[1]: if gu=FAIL then RETURN(FAIL): fi: print(`A Linear Recurrence Equation for The Diagonal Coefficients of The Power Series Of`): print(TOP/BOT^a): print(``): print(`By Shalosh B. Ekhad `): print(``): print(` Let `, A(n), `be the coefficient of`, mul(x[i]^n,i=1..nops(x)), `in the Maclaurin Expansion of the Formal Power series`, TOP/BOT^a ): print(`then `, A(n), `satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0): print(``): print(`Proof:`): print(`Thanks to Cauchy's Integral Formula, A(n) equals to `): print(` a constant (independent of n) times `): print(`the contour-integral with respect to the complex variables`, op(x)): print(`around any poly-circle around the origin, of the function`): print(F(n,op(x))=TOP/BOT^a/mul(x[i]^(n+1),i=1..nops(x))): print(`Let's cleverly constuct rational functions`): for i from 1 to nops(gu[2]) do print(R[i](n,op(x))=gu[2][i]): od: print(`With the motive that`): print(add(coeff(ope,N,i)*F(n+i,op(x)),i=0..degree(ope,N))= add(D[x[i]](R[i](n,op(x))/BOT^a/mul(x[j]^(n+1),j=1..nops(x))),i=1..nops(x))): print(``): print(`(Check!)`): print(`and the theorem follows upon contour-integrating with respect to`, op(x)): print(` QED! `): print(`This took`, time()-t0, `seconds . `): end: #DIAGpaperSR(TOP,BOT,a,x,n,A,n0,MaxC): Inputs polynomials in the list of variables #x, TOP, and BOT, and a numeric or symbolic number a, and a symbol m, and a set of parameters pars #(if you are lazy, you can make it {}) and letters F,A,R, (for formulating the proof) #outputs a mathematical article #authored by Shalosh B. Ekhad that conjectures but semi-rigorously proves a linear recurrence equation for the #diagonal sequence (called A(n)), i.e. the coefficient of x[1]^n*x[2]^n* in the formal power series #TOP/BOT^a. BOT must have a non-zero constant terms (for the formal power series to be well-defined) #For example, for Miklos Bona's problem, try: #DIAGpaperSR(1,(1-x)^2+(1-y)^2-1,1/2,[x,y],n,A,30,10); DIAGpaperSR:=proc(TOP,BOT,a,x,n,A,n0, MaxC) local gu,lu,i,j,mu,N,ope,t0: gu:=TOP/BOT^a: for i from 1 to nops(x) do gu:=taylor(gu,x[i],n0+1): gu:=add(coeff(gu,x[i],j)*x[i]^j,j=0..n0): od: lu:=[]: for i from 1 to n0 do mu:=gu: for j from 1 to nops(x) do mu:=coeff(mu,x[j],i): od: lu:=[op(lu),mu]: od: ope:=Findrec(lu,n,N,MaxC): if ope=0 or ope=FAIL then RETURN(FAIL): fi: print(`A Semi-Rigorously Proved Linear Recurrence Equation for The Diagonal Coefficients of The Power Series Of`): print(TOP/BOT^a): print(``): print(`By Shalosh B. Ekhad `): print(``): print(` Let `, A(n), `be the coefficient of`, mul(x[i]^n,i=1..nops(x)), `in the Maclaurin Expansion of the Formal Power series`, TOP/BOT^a ): print(`then `, A(n), `satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0): print(``): print(`Proof: By the Fundamental Theorem of Wilf-Zeilberger theory A(n) must satisfy SOME linear recurrence equation with polynomial coefficients.`): print(`Since this recurrence holds for the first`, n0, `terms, is is probably the right one!`): print(`This took`, time()-t0, `seconds . `): end: