###################################################################### ## GuessHolo2: Save this file as GuessHolo2. To use it, stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read GuessHolo2: # ## Then follow the instructions given there. # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg@math.rutgers.edu. # ###################################################################### with(SolveTools): print(`First Written: April 15, 2005: tested for Maple 10`): print(`Previous Version: Feb. 2, 2006. `): print(`Last version: June 5, 2006 [thanks to Moa Apagodu]`): print(): print(`This is GuessHolo2, A Maple package`): print(`accompanying Doron Zeilberger's article: `): print(`The Holonomic Ansatz: I. Foundations`): print(`It guesses recurrence equations for discrete variables of TWO variables`): 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(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(`For a list of the supporting functions type: ezra1();`): print(): ezra1:=proc() if args=NULL then print(` GuessHolo2: A Maple package for Handling Holnomic Functions`): print(`of TWO discrete variables. The supporting procedures are`): print(): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: AF, ApplyOper, AsyC`): print(`CheckCanSmnO ,CoeffP,CoeffPs,Comm, CommSeq, DiagOper,`): print(`Div1, EvalOpeF2 , findrec,Findrec, FindrecF, `): print(`GH2c, GH2dir, GH2Polc, GH2PolDir,GH2Mdir, `): print(`GH2Mdir1, GH2MPol, GH2MPolDir, GH2MPolPar`): print(` GR1, GR1Pol, GR1Rat, Kr, Mult, Mult1, SeqFromRec`): print(` Sp,Tp, Yafe1 `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` GuessHolo2: A Maple package for empirically guessing partial recurrence`): print(`equations satisfied by Discrete Functions of TWO Variables`): print(): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: `): print(` Asy, AsyC, BaW, BaWp, BaWg, CanOpe, DiagonalStory `): print(`GH2, GH2D,GH2DPol, GH2Pol`): print(`GH2PolPar,GuessDiag, GuessDiagPar,GuessDiagPol,GuesssDiagPolPar,`): print(` HypS, LatticePathStory, LaW `): print(`LaWp, ReciPolStory, Y3, Zinn`): print(): elif nargs=1 and args[1]=AF then print(`AF(F,n,k,m0,n0): Given a function F(n,k) finds`): print(`sum(F(n0+m0,k),k=0..m0)`): print(`For example, try: AF(binomial(n,k),n,k,3,0);`): elif nargs=1 and args[1]=ApplyOper then print(`ApplyOper(Seq1,ope,n,N,i0): Given a sequence Seq1, applies the operator `): print(`ope(n,N) at Seq1[i]`): elif nargs=1 and args[1]=Asy then print(`Asy(ope,n,N,K): the asymptotic expansion of solutions `): print(`to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator`): print(`up to the K's term. Note: ope has to be the MINIMAL operator!`): print(`For example, try: Asy((n^2+3)*N-(n^2+1),n,N,2);`): elif nargs=1 and args[1]=AsyC then print(`AsyC(ope,n,N,K,Ini,L): the asymptotic expansion `): print(`complete with the estimated constant, of solutions `): print(`to ope(n,N)f(n)=0, with initial conditions`): print(`Ini where ope(n,N) is a recurrence operator`): print(`up to the K's term. The constant is estimated using up to L terms`): print(`For example, try`): print(`AsyC(N^2-N-1,n,N,1,[1,1],100);`): elif nargs=1 and args[1]=BaW then print(`BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps ending at the point a`): print(`For example, try: BaW({[1,0],[0,1]},[1,1]);`): elif nargs=1 and args[1]=BaWp then print(`BaWp(Steps,a): The prob. of ballot-lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps, with associated probabilities, ending at the point a`): print(`For example, try: BaWp({[[1,0],1/2],[[0,1],1/2]},[1,1]);`): elif nargs=1 and args[1]=BaWg then print(`BaWg(Steps,a,g): The number of generalized ballot-lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps ending at the point a`): print(`These walks stay at a[i]-a[i+1]>=g`): print(`For example, try: BaWg({[1,0],[0,1]},[10,10],6);`): elif nargs=1 and args[1]=CanOpe then print(`CanOpe(Ope,m,n,M,N): Given a partial linear-recurrence operator with polynomial coefficients`): print(`divides it by an appropriate shiftp operaor to find an equivalent polynomial`): print(`in the form of a polynomial in (m,n,1/M,1/N)`): print(`For example, try: CanOpe(M*N-M-N,m,n,M,N);`): elif nargs=1 and args[1]=CanOpePos then print(`CanOpePos(Ope,m,n,M,N): Given a partial linear-recurrence operator with polynomial coefficients`): print(`converts it to `): print(`in the form of a polynomial in (m,n,M,N)`): print(`For example, try: CanOpePos(M*N-(1/m)*M-(1/n)*N,m,n,M,N);`): elif nargs=1 and args[1]=CanSmn then print(`CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representation`): print(`of a discrete function (let's call it F(m,n)), in terms`): print(`of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses`): print(`F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=abinomial(i+j,j),M*N-M-N,m,n,M,N,3,5);`): elif nargs=1 and args[1]=findrec then print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`): print(`the sequence f of degree DEGREE and order ORDER.`): print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`): print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `): print(`of maximum DEGREE+ORDER<=MaxC`): print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`): elif nargs=1 and args[1]=FindrecF then print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`): print(`poly coffs. .g. try FindrecF(i->i!,n,N);`): elif nargs=1 and args[1]=HypS then print(`HypS(F,n,k,m,M,Patience,s1): The Operator Oper(M,m,n) annihilating`): print(`sum(F(n+m,k),k=0..m)`): print(`For example, try: HypS(binomial(n,k),n,k,m,M,10,10);`): elif nargs=1 and args[1]=GH2 then print(`GH2(F,m,n,M,N,Patience,s1): Given a function of two variables (m,n) gives the holnomic`): print(`representation,[Ope1(m,n,M),Ope2(m,n,N)], where M and N`): print(`are the shift-operator in the m- and n- directions resp. `): print(`For example, try: GH2((i,j)->i!+j!,m,n,M,N,10,15);`): elif nargs=1 and args[1]=GH2dir then print(`GH2dir(F,m,n,M,N): Given a function of two variables (m,n) gives the holnomic`): print(`representation,[Ope1(m,n,M),Ope2(m,n,N)], where M and N`): print(`are the shift-operator in the m- and n- directions resp. `): print(`It uses the DIRECT Way`): print(`For example, try: GH2dir((i,j)->i!+j!,m,n,M,N);`): elif nargs=1 and args[1]=GH2c then print(`GH2c(F,m,n,M,N,Patience,s1): Given a function of two variables (m,n) gives the holnomic`): print(`representation,[Ope1(m,n,M),Ope2(m,n,N)], where M and N`): print(`are the shift-operator in the m- and n- directions resp. `): print(`In Canonical Form`): print(`For example, try: GH2c((i,j)->i!+j!,m,n,M,N,10,15);`): elif nargs=1 and args[1]=GH2D then print(`GH2D(F,m,n,M,N): Given a function of two variables (m,n) guesss a recurrence`): print(`in the diagonal direction (MN), where M and N are the shift-operators in the m-direction `): print(`and n-direction resp.`): print(`For example, try: GH2D((i,j)->i!+j!,m,n,M,N);`): elif nargs=1 and args[1]=GH2M then print(`GH2(F,m,n,M,Patience,s1):`): print(` Given a function of two variables (m,n) guesses a recurrence`): print(`in the m direction, where M is the shift-operator in the m-direction `): print(`Patience is the patience, and s1 is the starting point`): print(`For example, try: GH2M((i,j)->i!+j!,m,n,M,10,11);`): elif nargs=1 and args[1]=GH2Mdir then print(`GH2Mdir(F,m,n,M): Given a function of two variables (m,n) guesses a recurrence`): print(`in the m direction, where M is the shift-operator in the m-direction `): print(`For example, try: GH2Mdir((i,j)->i!+j!,m,n,M);`): elif nargs=1 and args[1]=GH2Mdir1 then print(`GH2Mdir1(F,m,n,M,deg,ORD), guesses an operator ope(m,n,M) of degree deg in (m,n)`): print(`and order ORD in M annihilating the function `): print(`for example, try`): print(`GH2Mdir1((i,j)->(i+j)!/i!/j!,m,n,M,1,1);`): elif nargs=1 and args[1]=GH2MPol then print(`GH2MPol(POL,x,y,m,n,M,Patience,s1): The recurrence in the m-direction of`): print(`the coeff. of x^m*y^n of 1/POL(x,y)`): print(`For example, try: GH2MPol(1-x-y,x,y,m,n,M,10,11);`): elif nargs=1 and args[1]=GH2MPolDir then print(`GH2MPolDir(POL,x,y,m,n,M,Patience): The recurrence in the m-direction of`): print(`the coeff. of x^m*y^n of 1/POL(x,y)`): print(`Doing it directly`): print(`For example, try: GH2MPolDir(1-x-y,x,y,m,n,M,10);`): elif nargs=1 and args[1]=GH2MPolPar then print(`GH2MPolPar(POL,x,y,m,n,M,par,Patience,s1): The recurrence in the m-direction of`): print(`the coeff. of x^m*y^n of 1/POL(x,y), where POL depends on par`): print(`For example, try: GH2MPolPar(1-x-y+a*x*y,x,y,m,n,M,a,10,11);`): elif nargs=1 and args[1]=GH2Pol then print(`GH2Pol(POL,x,y,m,n,M,N,Patience,s1): A holonomic representation`): print(`the coeff. of x^m*y^n of 1/POL(x,y).`): print(`Here M and N are the shift operators in the m- and n-directions`): print(`Patience is he pateice level (a positive integer >5)`): print(` s1 is the starting place of investigation`): print(`For example, try: GH2Pol(1-x-y,x,y,m,n,M,N,10,5);`): elif nargs=1 and args[1]=GH2PolDir then print(`GH2PolDir(POL,x,y,m,n,M,N): A holonomic representation`): print(`the coeff. of x^m*y^n of 1/POL(x,y).`): print(`Here M and N are the shift operators in the m- and n-directions`): print(`Patience is he pateice level (a positive integer >5)`): print(`For example, try: GH2PolDir(1-x-y,x,y,m,n,M,N);`): elif nargs=1 and args[1]=GH2PolPar then print(`GH2PolPar(POL,x,y,m,n,M,N,par,Patience,s1): A holonomic representation`): print(`the coeff. of x^m*y^n of 1/POL(x,y).`): print(`where POL depends on a parameter par.`): print(`Here M and N are the shift operators in the m- and n-directions`): print(`Patience is the patience, and s1 is the starting point`): print(`For example, try: GH2PolPar(1-x-a*y,x,y,m,n,M,N,a,20,11);`): elif nargs=1 and args[1]=GH2Polc then print(`GH2Polc(POL,x,y,m,n,M,N): A holonomic representation`): print(`the coeff. of x^m*y^n of 1/POL(x,y).`): print(`Here M and N are the shift operators in the m- and n-directions`): print(`In Canonical Form`): print(`For example, try: GH2Polc(1-x-y,x,y,m,n,M,N);`): elif nargs=1 and args[1]=GH2DPol then print(`GH2DPol(POL,x,y,m,n,M,N): The recurrence in the diagonal of`): print(`the coeff. of x^m*y^n of 1/POL(x,y).`): print(`Here M and N are the shift operators in the m- and n-directions`): print(`For example, try: GH2DPol(1-x-y,x,y,m,n,M,N);`): elif nargs=1 and args[1]=GR1 then print(`GR1(S1,x): guesses a polynomial from the data-set S1`): elif nargs=1 and args[1]=GR1Pol then print(`GR1Pol(S1,N,x): Guesses a poly in N rational in x for`): print(`the data set S1 of the form {[i,POL(N)]}:`): elif nargs=1 and args[1]=GR1Rat then print(`GR1Rat(S1,N,x): Guesses a rat. function in N rational in x for`): print(`the data set S1 of the form {[i,Rat(N)]}:`): elif nargs=1 and args[1]=GuessDiag then print(`GuessDiag(F,n,N,a): Given a function F of two discrete variables, guesses`): print(`a recurrence for the diag. F(n+a,n), where N is the shift in N`): print(`For example, try: GuessDiag((i,j)->i!+j!,n,N,0);`): elif nargs=1 and args[1]=GuessDiagPar then print(`GuessDiagPar(F,n,N,a,p): Given a function of two variables (m,n) `): print(`(featuring a parameter p), guesss a recurrence`): print(`in the diagonal (n+a,n)`): print(`For example, try: GuessDiagPar((i,j)->p*i!+j!,n,N,1,p);`): elif nargs=1 and args[1]=GuessDiagPol then print(`GuessDiagPol(P,x,y,n,N,a): Given a polynomial P of two variables x and y, `): print(`guesses a recurrence for the coeffs. of x^n*y^(n+a)`): print(`, where N is the shift in N.`): print(`For example, try: GuessDiagPol(1-x-y+2*x*y,x,y,n,N,0);`): elif nargs=1 and args[1]=GuessDiagPolPar then print(`GuessDiagPolPar(P,x,y,n,N,a): Given a polynomial P of two variables x and y, `): print(`That also depends (rationally) on a parameter p,`): print(`guesses a recurrence for the coeffs. of x^n*y^(n+a)`): print(`, where N is the shift in N.`): print(`For example, try: GuessDiagPolPar(1-x-y+p*x*y,x,y,n,N,0,p);`): elif nargs=1 and args[1]=Kr then print(`Kr(i,j)=Kronecker delta`): elif nargs=1 and args[1]=LatticePathStory then print(`LatticePathStory(Steps,Patience, s1): the story about Lattice Paths with Fundamental Steps`): print(`Steps: For example, try: LatticePathStory({[0,1],[1,0],[1,1]},10,15);`): elif nargs=1 and args[1]=LaW then print(`LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps ending at the point a`): print(`For example, try: LaW({[1,0],[0,1]},[1,1]);`): elif nargs=1 and args[1]=LaWp then print(`LaWp(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps, with associated probabilities ending at the point a`): print(`For example, try: LaWp({[[1,0],1/2],[[0,1],1/2]},[1,1]);`): elif nargs=1 and args[1]=Mult then print(`Mult(Ope1,Ope2,m,n,M,N): Ope1(m,n,M,N)*Ope2(m,n,M,N)`): elif nargs=1 and args[1]=Mult1 then print(`Mult1(ope1,ope2,n,N): ope1(n,N)*ope2(n,N)`): elif nargs=1 and args[1]=ReciPolStory then print(`ReciPolStory(P,x,y,Patience, s1): the story about the Maclaurin coefficients of the rational`): print(`function 1/P, where P is a polynomial in x,y`): print(`For example, try: ReciPolStory(1-x-y,x,y,10,15);`): elif nargs=1 and args[1]=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): print(`For example, try:`): print(`SeqFromRec(N-n-1,n,N,[1],10);`): elif nargs=1 and args[1]=Sp then print(`Sp(f,x,p): The inverse transformation to T_p(f(x)): for example, try`): print(` Sp(Tp(1/(1-x-y),[x,y],1/3),[x,y],1/3); `): elif nargs=1 and args[1]=Tp then print(`Tp(f,x,p): The transformation T_p(f(x)`) : print(`x[i]=p*x[i]/(1-(1-p)*x[i]), i=1..nops(x)`): print(`divided by mul(1-(1-p)*x[i],i=1..nops(x))`): print(`It preserves absolute `): print(`monotinicity for 0<=p<=1. Here f is a rational function in the list of variables `): print(`x, and p is beween 0 and 1 `): print(`for example, try`): print(` Tp(1/(1-x-y),[x,y],1/3); `): elif nargs=1 and args[1]=Y3 then print(`Y3(L): The number of 3-D Young Tableaux of shape L`): print(`For example, try: Y3([[2,2],[2,2]]);`): elif nargs=1 and args[1]=Yafe1 then print(`Yafe1(Ope,M,N): Makes Ope nice`): print(`For example, try: Yafe1((1/m)*M*N-1,M,N);`): elif nops([args])=1 and op(1,[args])=Zinn then print(` Zinn's method `): print(`Zinn(IntegerSequence): Given an increasing sequence a(n) of positive integers , expressed`): print(`in terms of a list, estimates the theta and mu such that `): print(` a(n) is asympt. to n^(theta)*mu^n. The output is `): print(` theta, mu `): fi: end: #####One variable rational-function geussing #Bdok1(DS1,P,x): checks that P in the single variable x indeed fits the data set DS1 Bdok1:=proc(DaS,P,x) local i: for i from 1 to nops(DaS) do if subs(x=DaS[i][1],denom(P))<>0 and normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then RETURN(false): fi: od: true: end: #GenRat1(x,dxt,dxb,a): a generic (monic top) #rational function of degree of top dxt #and degree of bottom dxb #with coeffs.parametrized by a[i] followed by the set #of coeffs. GenRat1:=proc(x,dxt,dxb,a) local i: (add(a[i]*x^i,i=0..dxt))/ add(a[i+dxt+1]*x^i,i=0..dxb) ,{seq(a[i],i=0..dxt+dxb+1)}: end: #GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2 #that fits the data in DS1 GR1d1d2:=proc(S1,x,d1,d2) local eq,var,P,a,var1,i: P:=GenRat1(x,d1,d2,a): var:=P[2]: P:=P[1]: if nops(S1)<=d1+d2+3 then # print(`Insufficient data`): RETURN(FAIL): fi: eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}: var1:=Linear(eq,var): if expand(subs(var1,numer(P)))=0 then RETURN(FAIL): fi: P:=subs(var1,P): if Bdok1(S1,P,x) then RETURN(factor(normal(P))): else RETURN(FAIL): fi: end: #GR1(S1,x): guesses a polynomial from the data-set S1 GR1:=proc(S1,x) local d1,d2,gu,L,i1: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: for L from 0 to nops(S1)-4 do for d1 from 0 to L do d2:=L-d1: gu:=GR1d1d2(S1,x,d1,d2): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: #GR1start(S1,x,D1,D2): guesses a polynomial from the data-set S1 #starting with degree D1/D2 GR1start:=proc(S1,x,D1,D2) local d1,d2,gu,L,i1: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: for d1 from D1 to nops(S1)-3-D2 do for d2 from D2 to nops(S1)-3-d1 do gu:=GR1d1d2(S1,x,d1,d2): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: #GR1Pol(S1,N,x): Guesses a poly in N rational in x for #the data set S1 of the form {[i,POL(N)]}: GR1Pol:=proc(S1,N,x) local gu,S1a,Deg,lu,i,i1,d1,d2,lDeg: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))): lDeg:=min(seq(ldegree(S1[i1][2],N),i1=1..nops(S1))): gu:=0: S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,lDeg)],i1=1..nops(S1))}: lu:=GR1(S1a,x): if lu=FAIL then RETURN(FAIL): else gu:=gu+factor(normal(lu)*N^lDeg): fi: d1:=degree(numer(lu),x): d2:=degree(denom(lu),x): for i from lDeg+1 to Deg do S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}: lu:=GR1start(S1a,x,d1,d2): if lu=FAIL then RETURN(FAIL): else gu:=gu+factor(normal(lu))*N^i: fi: od: gu: end: #GR1PolRat(S1,N,x,M): Guesses a poly in N rational in x for #the data set S1 of the form {[i,POL(N,Rat(M))]}: GR1PolRat:=proc(S1,N,x,M) local gu,S1a,Deg,lu,i,i1: option remember: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))): gu:=0: for i from 0 to Deg do S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}: lu:=GR1Rat(S1a,M,x): if lu=FAIL then RETURN(FAIL): else gu:=gu+factor(normal(lu))*N^i: fi: od: gu: end: ####End one-variable rational function guessing ###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 1 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: #findrecK(f,DEGREE,ORDER): Is there a unique recurrence if degree DEGREE and order ORDER #the sequence f of degree DEGREE and order ORDER #For example, try: findrecK([seq(i,i=1..10)],0,2); findrecK:=proc(f,DEGREE,ORDER) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,n,N,a: option remember: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: #if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then # RETURN(FAIL): #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(false): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)=1 then RETURN(true): else RETURN(false): fi: end: #end Findrec ####Begin Asymptotics #Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N) Aope1:=proc(ope,n,N) local gu: gu:=expand(numer(normal(ope))): gu:=coeff(gu,n,degree(gu,n)): gu:=expand(gu/coeff(gu,N,degree(gu,N))): factor(gu): end: C3:=proc(a,b,c) option remember: if [a,b,c]=[0,0,0] then RETURN(1): fi: if a<0 or b<0 or c<0 then RETURN(0): fi: C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1): end: #Aluf1(pit): is it positive dominant? Aluf1:=proc(pit) local aluf,shi,i,makom: shi:=evalf(abs(pit[1])): aluf:={evalf(evalc(pit[1]))}: makom:={1}: for i from 2 to nops(pit) do if evalf(abs(pit[i]))>evalf(shi) then shi:=abs(evalf(evalc(pit[i]))): aluf:={pit[i]}: makom:={i}: elif evalf(abs(pit[i]))=shi then aluf:=aluf union {evalf(evalc(pit[i]))}: makom:=makom union {i}: fi: od: aluf,shi,makom: end: #OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepAS1:=proc(ope1,n,N,alpha,f,S1) local x,f1,L,F,A,mu,i,A1: f1:=subs(n=1/x,f): L:=degree(f1,x): F:=f1+A*x^(L+S1): mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F) ,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): mu:=taylor(mu,x=0,L+11): mu:=simplify(mu): for i from L+1 to L+9 while coeff(mu,x,i)=0 do od: if i=L+10 then RETURN(FAIL): fi: mu:=coeff(mu,x,i): A1:=subs(Linear({mu},{A}),A): if A1=0 then RETURN(FAIL): fi: subs({x=1/n,A=A1},F): end: #OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepA:=proc(ope1,n,N,alpha,f) local S1: for S1 from 1 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od: OneStepAS1(ope1,n,N,alpha,f,S1): end: #Asy(ope,n,N,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the K's term Asy:=proc(ope,n,N,K) local gu,pit,lu,makom,alpha,mu,ope1,ku,i,f,x,ka: gu:=Aope1(ope,n,N): if degree(gu,N)=0 then print(`Not of exponential growth, case not implemented`): RETURN(FAIL): fi: if not type(expand(normal(ope/gu)),`+`) then pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: RETURN(pit[makom[1]]^n*expand((nops(makom)-1)!*binomial(nops(makom)-1+n,nops(makom)-1))): fi: pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Not only that but Dominant roots are complex`): RETURN(FAIL): fi: fi: lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Dominant root is complex`): RETURN(FAIL): fi: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,nops(makom)+1),x,nops(makom))): ka:=simplify([solve(ku,alpha)]): alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN(lu^n*n^alpha): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: lu^n*n^alpha*f: end: ###EndAsymptotics #####Coeffs. of reciprocal of a polynomial #PolToSteps(POL,xList): Inputs a polynomial POL in #the set of variables xList and returns the #Weighted Steps for the scheme followed by the reciprocal #of the constant term PolToSteps:=proc(POL,xList) local k,i,POL1,mu,Steps,ku,ku1,i1,pu: option remember: POL1:=expand(POL): k:=nops(xList): if not type(POL,`+`) then print(`POL must have at least two terms`): RETURN(FAIL): fi: pu:=POL1: for i from 1 to nops(xList) do pu:=coeff(pu,xList[i],0): od: if pu=0 then print(`The polynomial should have a non-zero constant term`): RETURN(FAIL): fi: POL1:=expand(POL1/pu): Steps:={}: for i from 1 to nops(POL1) do mu:=op(i,POL1): ku:=[seq(degree(mu,xList[i1]),i1=1..nops(xList))]: ku1:=normal(mu/convert([seq(xList[i1]^ku[i1],i1=1..nops(xList))],`*`)): Steps:=Steps union {[ku,-ku1]}: od: Steps:=Steps minus {[[0$k],-1]}: Steps,1/pu: end: #Coeff1(Steps,a): #The number of weighter-lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps ending at the point a #For example, try: Coeff1({[[1,0],1],[[0,1],1]},[1,1]); Coeff1:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: Prev:= {seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}: expand(add(Prev[i][2]*Coeff1(Steps,Prev[i][1]),i=1..nops(Prev))): end: #Coeff1S(Steps,a): #The number of weighted-lattice walks in the k-dim lattice with pos. steps #in the set of symmetric steps Steps ending at the point a #For example, try: Coeff1({[[1,0],1],[[0,1],1]},[1,1]); Coeff1S:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: Prev:= {seq([[seq(a[j]-Steps[i][ 1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}: expand(add(Prev[i][2]*Coeff1S(Steps,sort(Prev[i][1])),i=1..nops(Prev))): end: #CoeffP(POL,xList,aList): The coeff. of xList^aList in POL #(a polynomial in the variables of xList) CoeffP:=proc(POL,xList,aList) local Steps,coe: Steps:=PolToSteps(POL,xList): coe:=Steps[2]: Steps:=Steps[1]: coe*Coeff1(Steps,aList): end: #CoeffPs(POL,xList,aList): The coeff. of xList^aList in 1/POL #(a Symmetric polynomial in the variables of xList) CoeffPs:=proc(POL,xList,aList) local Steps,coe: Steps:=PolToSteps(POL,xList): coe:=Steps[2]: Steps:=Steps[1]: coe*Coeff1S(Steps,aList): end: ##### #LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps ending at the point a #For example, try: LaW({[1,0],[0,1]},[1,1]); LaW:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: add(LaW(Steps,Prev[i]),i=1..nops(Prev)): end: #LaWp(Steps,a): The probability of lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps with attached probabilities #to get to the point a #For example, try: LaWp({[[1,0],1/2],[[0,1],1/2]},[1,1]); LaWp:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: Prev:={seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}: add(Prev[i][2]*LaWp(Steps,Prev[i][1]),i=1..nops(Prev)): end: #BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps ending at the point a #For example, try: BaW({[1,0],[0,1]},[1,1]); BaW:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then RETURN(0): fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: add(BaW(Steps,Prev[i]),i=1..nops(Prev)): end: #BaWp(Steps,a): The probability of ballot-lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps, with assoc. prob. ending at the point a #For example, try: BaWp({[[1,0],1/2],[[0,1],1/2]},[1,1]); BaWp:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then RETURN(0): fi: Prev:={seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}: add(Prev[i][2]*BaWp(Steps,Prev[i][1]),i=1..nops(Prev)): end: #BaWg(Steps,a,g): The number of generalized ballot-lattice walks #(i.e. that stay in a[i]-a[i+1]>=-g #in the k-dim lattice with pos. steps #in the set of steps Steps ending at the point a #For example, try: BaWg({[1,0],[0,1]},[1,1],2); BaWg:=proc(Steps,a,g) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: if min(seq(a[i]-a[i+1]+g,i=1..nops(a)-1))<0 then RETURN(0): fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: add(BaWg(Steps,Prev[i],g),i=1..nops(Prev)): 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: #FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with #poly coffs. #e.g. try FindrecF(i->i!,n,N); FindrecF:=proc(f,n,N) local DEGREE, ORDER,ope,L,i,t0: option remember: t0:=time(): for L from 0 do for ORDER from 0 to L do DEGREE:=L-ORDER: ope:=findrec([seq(f(i),i=1..(2+DEGREE)*(1+ORDER)+4)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #GH2MPol(POL,x,y,m,n,M,Patience,s1): The recurrence in the m-direction of #the coeff. of x^m*y^n of 1/POL(x,y) #For example, try: GH2MPol(1-x-y,x,y,m,n,M,10); GH2MPol:=proc(P,x,y,m,n,M,Patience,s1) local i,j: if expand(subs({x=y,y=x},P)-P)=0 then GH2M((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,Patience,s1): else GH2M((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,Patience,s1): fi: end: GH2:=proc(F,m,n,M,N,Patience,s1) local i,j: GH2M(F,m,n,M,Patience,s1),GH2M((i,j)->F(j,i),n,m,N,Patience,s1):end: GH2dir:=proc(F,m,n,M,N,MaxC) local i,j: GH2Mdir(F,m,n,M),GH2Mdir((i,j)->F(j,i),n,m,N):end: GH2PolDir:=proc(P,x,y,m,n,M,N) local ope: if expand(subs({x=y,y=x},P)-P)=0 then ope:=GH2Mdir((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M): RETURN(ope,subs({m=n,n=m,M=N},ope)): else GH2dir((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N): fi: end: #GH2c(F,m,n,M,N,Patience,s1): Canonical Form of GH2 GH2c:=proc(F,m,n,M,N,Patience,s1) local lu,lu1,lu2: lu:=GH2(F,m,n,M,N,Patience,s1): lu1:=numer(normal(lu[1])):lu2:=numer(normal(lu[2])): Yafe1(lu1,M,N),Yafe1(lu2,M,N):end: GH2Polc:=proc(P,x,y,m,n,M,N)local lu: lu:=GH2c((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N): Yafe1(lu[1],M,N),Yafe1(lu[2],M,N):end: GH2DPol:=proc(P,x,y,m,n,M,N):GH2D((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):end: #EvalOpeF2(F,ope,m,n,M,N,m0,n0): evaluates ope(m,n,M,N)F(m0,n0) #For example try: #EvalOpeF2((i,j)->binomial(i+j,j),M*N-M-N,m,n,M,N,3,5); EvalOpeF2:=proc(F,ope,m,n,M,N,m0,n0) local gu,i,j: gu:=0: for i from 0 to degree(ope,M) do for j from 0 to degree(ope,N) do gu:=gu+subs({m=m0,n=n0},coeff(coeff(ope,M,i),N,j))*F(m0+i,n0+j): od: od: gu: end: #CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representation #of a discrete function (let's call it F(m,n)), in terms #of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses #F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a1 then print(OpeM, `should be monic in `, M): RETURN(FAIL): fi: if coeff(OpeN,N,Ordn)<>1 then print(OpeN, `should be monic in `, N): RETURN(FAIL): fi: if i0=Ordn if i0=Ordm and j0>=Ordn gu:=subs(m=m+1,CanSmn(OpeM,OpeN,m,n,M,N,i0-1,j0,F)): for j1 from 0 to Ordn-1 do gu:=expand(subs(F(m+Ordm,n+j1)=CanSmn(OpeM,OpeN,m,n,M,N,Ordm,j1,F),gu)): od: gu1:=0: for i1 from 0 to Ordm-1 do for j1 from 0 to Ordn-1 do gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*F(m+i1,n+j1): od: od: RETURN(gu1): end: #CanSmnO(OpeM,OpeN,m,n,M,N,i0,j0): #CanSmnO(M-1,N-1,m,n,M,N,1,1); CanSmnO:=proc(OpeM,OpeN,m,n,M,N,i0,j0) local F,i1,j1,gu,gu1,Ordm,Ordn: gu:=CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Ordm:=degree(OpeM,M): Ordn:=degree(OpeN,N): gu1:=0: for i1 from 0 to Ordm-1 do for j1 from 0 to Ordn-1 do gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1: od: od: gu1:=M^i0*N^j0-gu1: RETURN(gu1): end: #CheckCanSmnO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P #For example: #CheckCanSmnO(1-x-y-x*y,x,y,3,4); CheckCanSmnO:=proc(P,x,y,i0,j0) local lu,m,n,M,N,Ope,i,j,m0,n0: lu:=GH2Pol(P,x,y,m,n,M,N,5): Ope:=CanSmnO(lu[1],lu[2],m,n,M,N,i0,j0): evalb( {seq(seq(EvalOpeF2((i,j)->CoeffP(P,[x,y],[i,j]),Ope,m,n,M,N,m0,n0),m0=0..40),n0=0..40)} ={0}) : end: #CanSmnO1(OpeM,OpeN,m,n,M,N,i0,j0): #CanSmnO1(M-1,N-1,m,n,M,N,1,1); CanSmnO1:=proc(OpeM,OpeN,m,n,M,N,i0,j0) local F,i1,j1,gu,gu1,Ordm,Ordn: gu:=CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Ordm:=degree(OpeM,M): Ordn:=degree(OpeN,N): gu1:=0: for i1 from 0 to Ordm-1 do for j1 from 0 to Ordn-1 do gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1: od: od: gu1:=gu1: RETURN(gu1): end: #CanSnm(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representation #of a discrete function (let's call it F(m,n)), in terms #of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses #F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a1 then print(OpeM, `should be monic in `, M): RETURN(FAIL): fi: if coeff(OpeN,N,Ordn)<>1 then print(OpeN, `should be monic in `, N): RETURN(FAIL): fi: if i0=Ordn if i0=Ordm and j0>=Ordn gu:=subs(n=n+1,CanSmn(OpeM,OpeN,m,n,M,N,i0,j0-1,F)): for i1 from 0 to Ordm-1 do gu:=expand(subs(F(m+i1,n+Ordn)=CanSnm(OpeM,OpeN,m,n,M,N,i1,Ordn,F),gu)): od: gu1:=0: for i1 from 0 to Ordm-1 do for j1 from 0 to Ordn-1 do gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*F(m+i1,n+j1): od: od: RETURN(gu1): end: #CanSnmO(OpeM,OpeN,m,n,M,N,i0,j0): #CanSnmO(M-1,N-1,m,n,M,N,1,1); CanSnmO:=proc(OpeM,OpeN,m,n,M,N,i0,j0) local F,i1,j1,gu,gu1,Ordm,Ordn: gu:=CanSnm(OpeM,OpeN,m,n,M,N,i0,j0,F): Ordm:=degree(OpeM,M): Ordn:=degree(OpeN,N): gu1:=0: for i1 from 0 to Ordm-1 do for j1 from 0 to Ordn-1 do gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1: od: od: gu1:=M^i0*N^j0-gu1: RETURN(gu1): end: #CheckCanSnmO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P #For example: #CheckCanSnmO(1-x-y-x*y,x,y,3,4); CheckCanSnmO:=proc(P,x,y,i0,j0) local lu,m,n,M,N,Ope,i,j,m0,n0: lu:=GH2Pol(P,x,y,m,n,M,N,5): Ope:=CanSnmO(lu[1],lu[2],m,n,M,N,i0,j0): evalb( {seq(seq(EvalOpeF2((i,j)->CoeffP(P,[x,y],[i,j]),Ope,m,n,M,N,m0,n0),m0=0..40),n0=0..40)} ={0}) : end: #Yafe1(Ope,M,N): Makes Ope nice Yafe1:=proc(Ope,M,N) local d1,d2,Ope1,i1,i2,gu,Ld1,Ld2: #Ope1:=numer(normal(Ope)): Ope1:=Ope: d1:=degree(Ope1,M): Ld1:=ldegree(Ope1,M): d2:=degree(Ope1,N): Ld2:=ldegree(Ope1,N): gu:=0: for i1 from Ld1 to d1 do for i2 from Ld2 to d2 do gu:=gu+factor(coeff(coeff(Ope1,M,i1),N,i2))*M^i1*N^i2: od: od: gu: end: #Mult(Ope1,Ope2,m,n,M,N): Ope1(m,n,M,N)*Ope2(m,n,M,N) Mult:=proc(Ope1,Ope2,m,n,M,N) local d1,d2,i1,i2,gu,Ld1,Ld2: d1:=degree(Ope1,M): Ld1:=ldegree(Ope1,M): d2:=degree(Ope1,N): Ld2:=ldegree(Ope1,N): gu:=0: for i1 from Ld1 to d1 do for i2 from Ld2 to d2 do gu:=expand(gu+factor(coeff(coeff(Ope1,M,i1),N,i2))*M^i1*N^i2 *subs({m=m+i1,n=n+i2},Ope2) ): od: od: Yafe1(gu,M,N): end: #Comm(Ope1,Ope2,m,n,M,N): The commutator of Ope1, Ope2 #i.e. Ope1*Ope2-Ope2*Ope1 Comm:=proc(Ope1,Ope2,m,n,M,N) Yafe1(Mult(Ope1,Ope2,m,n,M,N)-Mult(Ope2,Ope1,m,n,M,N),M,N): end: #CommSeq(OPE0,Ope1,m,n,M,N): Given a partial recurrence operator #with constant coefficients, OPE0, and a partial recurrence #operator with polynomial coefficients, finds the #sequence of commutators with OPE0: [Ope1,[OPE0,Ope1],[OPE0,%], ...] #and checks that the last one is a multiple of OPE0. Otherwise #it returns FAIL CommSeq:=proc(OPE0,Ope1,m,n,M,N) local gu,lu,gu1: lu:=Ope1: gu:=[lu]: while lu<>0 do lu:=Comm(OPE0,lu,m,n,M,N): gu:=[op(gu),lu]: od: gu:=[op(1..nops(gu)-1,gu)]: gu:=[op(1..nops(gu)-1,gu),factor(gu[nops(gu)])]: gu1:=denom(ormal(gu[nops(gu)]/OPE0)); if gu1<>1 then FAIL: else gu: fi: end: Diag1:=proc(OpeM,OpeN,m,n,M,N,k0) local gu,a,var,eq,eqTry,i1,j1,OpeD,var1T,i,var1: OpeD:=M^k0*N^k0: var:={}: gu:=CanSmnO1(OpeM,OpeN,m,n,M,N,k0,k0): for i from 0 to k0-1 do gu:=gu+a[i]*CanSmnO1(OpeM,OpeN,m,n,M,N,i,i): var:= var union {a[i]}: OpeD:=OpeD+a[i]*M^i*N^i: od: gu:=expand(gu): eq:={seq(seq(coeff(coeff(gu,M,i1),N,j1),j1=0..degree(gu,N)), i1=0..degree(gu,M))}: eqTry:=subs({m=11/5,n=17/19},eq): var1T:=Linear(eqTry,var): if var1T=NULL then RETURN(FAIL): fi: var1:=Linear(eq,var): if var1=NULL then RETURN(FAIL): else RETURN(subs(var1,OpeD)): fi: end: #DiagOper(OpeM,OpeN,m,n,M,N):The operator in the MN direction #deduced from OpeM(M,m,n) and OpeN(N,m,n) #For example, try: DiagOper(M-1,N-1,m,n,M,N); DiagOper:=proc(OpeM,OpeN,m,n,M,N) local k0,gu: for k0 from 1 do gu:=Diag1(OpeM,OpeN,m,n,M,N,k0): if gu<>FAIL then RETURN(gu): fi: od: end: #GH2MPolPar(POL,x,y,m,n,M,par,Patience,s1): The recurrence in the m-direction of #the coeff. of x^m*y^n of 1/POL(x,y) #where POL also depends (in a rational way) #on P #For example, try: GH2MPolPar(1-x-y-a*x*y,x,y,m,n,M,a,10,11); GH2MPolPar:=proc(P,x,y,m,n,M,par,Patience,s1) local S1,i,Ope,Ope1,i1,j1,m0,n0,N: S1:={}: for i from 2 to 5 do S1:=S1 union {[i,GH2MPol(subs(par=i,P),x,y,m,n,M,Patience,s1)]}: od: for i from 6 do Ope:= GR1Pol(S1,M,par): if Ope<>FAIL then Ope1:=numer(normal(Ope)): if {seq(seq(normal(EvalOpeF2((i1,j1)->CoeffP(P,[x,y],[i1,j1]), Ope1,m,n,M,N,m0,n0)),m0=0..10),n0=0..10)}<>{0} then print( {seq(seq(normal(EvalOpeF2((i1,j1)->CoeffP(P,[x,y],[i1,j1]), Ope1,m,n,M,N,m0,n0)),m0=0..10),n0=0..10)}): print(Ope, ` failed `): else RETURN(Ope): fi: fi: S1:=S1 union {[i,GH2MPol(subs(par=i,P),x,y,m,n,M,Patience,s1)]}: od: end: #GH2PolPar(P,x,y,m,n,M,N,par,Patience,s1): Like GH2Pol, but P depends #on a parameter par GH2PolPar:=proc(P,x,y,m,n,M,N,par,Patience,s1) local ope: if expand(subs({x=y,y=x},P)-P)=0 then ope:=GH2MPolPar(P,x,y,m,n,M,par,Patience,s1): RETURN(ope,subs({m=n,n=m,M=N},ope)): else RETURN(GH2MPolPar(P,x,y,m,n,M,par,Patience,s1), subs({m=n,n=m,M=N},GH2MPolPar(subs({x=y,y=x},P),x,y,m,n,M,par,Patience,s1))) : fi: end: #GuessDiag(F,n,N,a): Given a function F of two discrete variables, guesses #a recurrence for the diag. F(n+a,n), where N is the shift in N #For example, try: GuessDiag((i,j)->i!+j!,n,N,0); GuessDiag:=proc(F,n,N,a) local i: FindrecF(i->F(i+a,i),n,N): end: #GuessDiagPol(P,x,y,n,N,a): Given a polynomial of two variables x and y, #guesses a recurrence for the coeffs. of x^n*y^(n+a) #, where N is the shift in N #For example, try: GuessDiagPol(1-x-y+2*x*y,x,y,n,N,0); GuessDiagPol:=proc(P,x,y,n,N,a) local i,j: GuessDiag((i,j)->CoeffP(P,[x,y],[i,j]),n,N,a): end: #GuessDiagPar(F,n,N,a,p): Given a function of two variables (m,n) #(featuring a parameter p), guesss a recurrence #in the diagonal (n+a,n) #For example, try: GuessDiagPar((i,j)->p*i!+j!,n,N,1,p); GuessDiagPar:=proc(F,n,N,a,p) local i,ope,h,ope1,ope2,ope3,DEG,ORD,gu1,gu2,gu3,S1, Ope,j: ope[1]:=FindrecF(h->subs(p=1,F(h+a,h)),n,N): ope[2]:=FindrecF(h->subs(p=2,F(h+a,h)),n,N): ope[3]:=FindrecF(h->subs(p=3,F(h+a,h)),n,N): ope1:=numer(normal(ope[1])): ope2:=numer(normal(ope[2])): ope3:=numer(normal(ope[3])): for i from 4 do ope1:=numer(normal(ope[i-1])): ope2:=numer(normal(ope[i-2])): ope3:=numer(normal(ope[i-3])): if nops({degree(ope1,n),degree(ope2,n),degree(ope2,n)})=1 and nops({degree(ope1,N),degree(ope2,N),degree(ope2,N)})=1 then DEG:=degree(ope3,n): ORD:=degree(ope3,N): gu1:=[seq(subs(p=i-1,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu2:=[seq(subs(p=i-2,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu3:=[seq(subs(p=i-3,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and findrecK(gu3,DEG,ORD) then break: fi: fi: ope[i]:=FindrecF(h->subs(p=i,F(h+a,h)),n,N): od: S1:={}: for j from i-2 to i+6 do gu1:=[seq(subs(p=j,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) then ope1:=findrec(gu1,DEG,ORD,n,N): S1:=S1 union {[j,ope1]}: fi: od: for j from i+7 do Ope:= GR1Pol(S1,N,p): if Ope<>FAIL then RETURN(Ope): fi: gu1:=[seq(subs(p=j,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) then ope1:=findrec(gu1,DEG,ORD,n,N): S1:=S1 union {[j,ope1]}: fi: od: end: #GuessDiagPolPar(P,x,y,n,N,a,p): Given a polynomial of two variables x and y, #also depending (rationally) on a parameter p #guesses a recurrence for the coeffs. of x^n*y^(n+a) #, where N is the shift in N #For example, try: GuessDiagPolPar(1-x-y+p*x*y,x,y,n,N,0,p); GuessDiagPolPar:=proc(P,x,y,n,N,a,p): GuessDiagPar((i,j)->CoeffP(P,[x,y],[i,j]),n,N,a,p): end: #GR1Rat(S1,N,x): Guesses a rational function in N rational in x for #the data set S1 of the form {[i,RAT(N)]}: GR1Rat:=proc(S1,N,x) local Deg1,Deg2,i,i1,S1Top,S1Bot,R1,c,Top,Bot: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: Deg1:=max(seq(degree(numer(S1[i1][2]),N),i1=1..nops(S1))): Deg2:=max(seq(degree(denom(S1[i1][2]),N),i1=1..nops(S1))): S1Top:={}: S1Bot:={}: for i from 1 to nops(S1) do R1:=S1[i][2]: if degree(numer(R1),N)=Deg1 and degree(denom(R1),N)=Deg2 then c:=coeff(numer(R1),N,Deg1): S1Top:=S1Top union {[S1[i][1],numer(R1)/c]}: S1Bot:=S1Bot union {[S1[i][1],denom(R1)/c]}: fi: od: Top:=GR1Pol(S1Top,N,x): Bot:=GR1Pol(S1Bot,N,x): if Top<>FAIL and Bot<>FAIL then RETURN(Top/Bot): else RETURN(FAIL): fi: end: #####Direct Way for GH2M: GH2Mdirect #GenPol2(x,y,d,a): a generic Polynomial #of dgree d in (x,y) with coeffs. indexed by a #followed by the coefficient GenPol2(x,y,1,a) #For example, try GenPol2:=proc(x,y,d,a) local i,k: convert([seq(seq(a[i,k-i]*x^i*y^(k-i),i=0..k),k=0..d)],`+`), {seq(seq(a[i,k-i],i=0..k),k=0..d)}: end: #GenOper2(m,n,M,deg,ORD,a): a generic operator of degree deg in (m,n) #and order ORD in M GenOper2:=proc(m,n,M,deg,ORD,a) local i,var,gu,lu: gu:=0: var:={}: for i from 0 to ORD do lu:=GenPol2(m,n,deg,a[i]): gu:=gu+lu[1]*M^i: var:=var union lu[2]: od: gu,var: end: #GH2Mdir1(F,m,n,M,deg,ORD), guesses an operator ope(m,n,M) of degree deg in (m,n) #and order ORD in M annihilating the function F #for example, try #GH2Mdir1((i,j)->(i+j!,m,n,M,2,2); GH2Mdir1:=proc(F,m,n,M,deg,ORD) local ope,var,eq,a,m0,n0,k,N,var1,ind1,i: ope:=GenOper2(m,n,M,deg,ORD,a): var:=ope[2]: ope:=ope[1]: eq:={}: for k from 5 while nops(eq)<=nops(var)+5 do for m0 from 0 to k do n0:=k-m0: eq:=eq union {EvalOpeF2(F,ope,m,n,M,N,m0+2,n0+4)}: od: od: var1:=Linear(eq,var): ind1:={}: for i from 1 to nops(var1) do if op(1,var1[i])=op(2,var1[i]) then ind1:=ind1 union {op(1,var1[i])}: fi: od: if nops(ind1)=0 then RETURN(FAIL): fi: if nops(ind1)>1 then RETURN(FAIL): fi: ope:=subs(ind1[1]=1,subs(var1,ope)): if {seq(seq(EvalOpeF2(F,ope,m,n,M,N,m0,n0),m0=10..20),n0=10..20)}<>{0} then RETURN(FAIL): else RETURN(Yafe1(ope,M,N)): fi: end: #GH2Mdir(F,m,n,M), guesses an operator ope(m,n,M) of degree deg in (m,n) #and order ORD in M annihilating the function F #for example, try #GH2Mdir((i,j)->(i+j!,m,n,M); GH2Mdir:=proc(F,m,n,M) local ope,DEG,ORD,k,LT,i,ope1,ope2,ope3,ope4,j: ope1:=numer(normal(FindrecF(i->F(i,20),n,N))): ope2:=numer(normal(FindrecF(i->F(i,21),n,N))): ope3:=numer(normal(FindrecF(i->F(i,22),n,N))): for j from 23 do if nops({degree(ope1,N),degree(ope2,N),degree(ope3,N)})=1 and nops({degree(ope1,n),degree(ope2,n),degree(ope3,n)})=1 then DEG:=degree(ope3,n): ORD:=degree(ope3,N): break: fi: ope4:=numer(normal(FindrecF(i->F(i,j),n,N))): ope1:=ope2: ope2:=ope3: ope3:=ope4: od: ope:=GH2Mdir1(F,m,n,M,DEG,ORD): if ope<>FAIL then LT:=coeff(ope,M,degree(ope,M)): ope:=add(normal(coeff(ope,M,i)/LT)*M^i,i=0..degree(ope,M)): RETURN(ope): fi: FAIL: end: #GH2MPolDir(POL,x,y,m,n,M,MaxC): The recurrence in the m-direction of #the coeff. of x^m*y^n of 1/POL(x,y) the direct way #For example, try: GH2MPolDir(1-x-y,x,y,m,n,M,10); GH2MPolDir:=proc(P,x,y,m,n,M,MaxC) local i,j: if expand(subs({x=y,y=x},P)-P)=0 then GH2Mdir((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,MaxC): else GH2Mdir((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,MaxC): fi: end: #GH2D(F,m,n,M,N): Given a function of two variables (m,n) guesss a recurrence #in the diagonal direction (MN), where M and N are the shift-operator in the m-direction #and n-directions #For example, try: GH2D((i,j)->i!+j!,m,n,M,N); GH2D:=proc(F,m,n,M,N) local i,ope,h,ope1,ope2,ope3,DEG,ORD,gu1,gu2,gu3,S1, Ope,j,K,x,Ope1,m0,n0,i1: ope[1]:=FindrecF(h->F(h+1,h),n,K): ope[2]:=FindrecF(h->F(h+2,h),n,K): ope[3]:=FindrecF(h->F(h+3,h),n,K): ope1:=numer(normal(ope[1])): ope2:=numer(normal(ope[2])): ope3:=numer(normal(ope[3])): for i from 4 do ope1:=numer(normal(ope[i-1])): ope2:=numer(normal(ope[i-2])): ope3:=numer(normal(ope[i-3])): if nops({degree(ope1,n),degree(ope2,n),degree(ope2,n)})=1 and nops({degree(ope1,K),degree(ope2,K),degree(ope2,K)})=1 then DEG:=degree(ope3,n): ORD:=degree(ope3,K): gu1:=[seq(F(h+i-1,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu2:=[seq(F(h+i-2,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu3:=[seq(F(h+i-3,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and findrecK(gu3,DEG,ORD) then break: fi: fi: ope[i]:=FindrecF(h->F(h+i,h),n,K): od: S1:={}: for j from i-2 to i+6 do gu1:=[seq(F(h+j,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) then ope1:=findrec(gu1,DEG,ORD,n,K): S1:=S1 union {[j,ope1]}: fi: od: for j from i+7 do Ope:= GR1PolRat(S1,K,x,n): if Ope<>FAIL then Ope:=subs({x=m-n},Ope): Ope:=add(factor(normal(coeff(Ope,K,i1)))*K^i1,i1=0..degree(Ope,K)): Ope:=subs({K=M*N},Ope): Ope1:=numer(normal(Ope)): if {seq(seq(normal(EvalOpeF2(F,Ope1,m,n,M,N,m0,n0)),m0=0..40),n0=0..40)}<>{0} then print(Ope, ` failed `): else RETURN(Ope): fi: fi: gu1:=[seq(F(h,j),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) then ope1:=findrec(gu1,DEG,ORD,n,K): S1:=S1 union {[j,ope1]}: fi: od: end: #GH2M(F,m,n,M,Patience,s1): #Given a function of two variables (m,n) guesses a recurrence #in the m direction, where M is the shift-operator in the m-direction #It does it the fast way (using GR1Rat instead of GR1) #s1 is where it starts exploring #For example, try: GH2M((i,j)->i!+j!,m,n,M,10); GH2M:=proc(F,m,n,M,Patience,s1) local i,ope,h,ope1,ope2,ope3,ope4,DEG,ORD,gu1,gu2,gu3,gu4,S1,Ope,j, m0,n0,N,Ope1: ope[1]:=FindrecF(h->F(h,s1+1),m,M): ope[2]:=FindrecF(h->F(h,s1+2),m,M): ope[3]:=FindrecF(h->F(h,s1+3),m,M): ope[4]:=FindrecF(h->F(h,s1+4),m,M): ope1:=numer(normal(ope[1])): ope2:=numer(normal(ope[2])): ope3:=numer(normal(ope[3])): ope4:=numer(normal(ope[4])): for i from 5 to Patience do ope1:=numer(normal(ope[i-1])): ope2:=numer(normal(ope[i-2])): ope3:=numer(normal(ope[i-3])): ope4:=numer(normal(ope[i-4])): if nops({degree(ope1,m),degree(ope2,m),degree(ope3,m),degree(ope4,m)})=1 and nops({degree(ope1,M),degree(ope2,M),degree(ope3,M),degree(ope4,M)})=1 then DEG:=degree(ope4,m): ORD:=degree(ope4,M): gu1:=[seq(F(h,s1+i-1),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu2:=[seq(F(h,s1+i-2),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu3:=[seq(F(h,s1+i-3),h=1..(1+DEG)*(1+ORD)+5+ORD)]: gu4:=[seq(F(h,s1+i-4),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and findrecK(gu3,DEG,ORD) and findrecK(gu4,DEG,ORD) then break: fi: fi: ope[i]:=FindrecF(h->F(h,i+s1),m,M): od: if i-1=Patience then print(`Couldn't do it with patience`, Patience, `and Starting place`, s1): RETURN(FAIL): fi: S1:={}: for j from i-3 to i+6 do gu1:=[seq(F(h,j+s1),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) then ope1:=findrec(gu1,DEG,ORD,m,M): S1:=S1 union {[j+s1,ope1]}: fi: od: if nops( {seq(S1[i][2],i=1..nops(S1))})=1 then RETURN(S1[1][2]): fi: for j from i+7 do Ope:= GR1PolRat(S1,M,n,m): if Ope<>FAIL then Ope1:=numer(normal(Ope)): if {seq(seq(normal(EvalOpeF2(F,Ope1,m,n,M,N,m0,n0)),m0=0..40),n0=0..40)}<>{0} then print(Ope, ` failed `): RETURN(FAIL): else RETURN(Ope): fi: fi: gu1:=[seq(F(h,j+s1),h=1..(1+DEG)*(1+ORD)+5+ORD)]: if findrecK(gu1,DEG,ORD) then ope1:=findrec(gu1,DEG,ORD,m,M): S1:=S1 union {[j+s1,ope1]}: fi: od: end: GH2Pol:=proc(P,x,y,m,n,M,N,Patience,s1) local ope: if expand(subs({x=y,y=x},P)-P)=0 then ope:=GH2M((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,Patience,s1): RETURN(ope,subs({m=n,n=m,M=N},ope)): else GH2M((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,Patience,s1), GH2M((i,j)->CoeffP(P,[x,y],[j,i]),m,n,N,Patience,s1) fi: end: #AF(F,n,k,m0,n0): Given a function F(n,k) finds #sum(F(n0+m0,k),k=0..m0) AF:=proc(F,n,k,m0,n0) option remember: if m0<0 then 0: else AF(F,n,k,m0-1,n0+1)+eval(subs({k=m0,n=m0+n0},F)): fi: end: #HypS(F,n,k,m,M,Patience,s1): The Operator Oper(M,m,n) annihilating #sum(F(n+m,k),k=0..m) #For example, try: HypS(binomial(n,k),n,k,m,M,10,10); HypS:=proc(F,n,k,m,M,Patience,s1): GH2M((i,j)->AF(F,n,k,i,j),m,n,M,Patience,s1): end: Zinn:=proc(resh) local s1,s2,n: n:=nops(resh)-2: s1:=sn(resh,n): s2:=sn(resh,n-1): evalf(2*(s1+s2)/(s1-s2)^2), evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): end: sn:=proc(resh,n): -1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2): #evalf(%): end: #CanOpe(Ope,m,n,M,N): Given a partial linear-recurrence operator with polynomial coefficients #divides it by an appropriate shiftp operaor to find an equivalent polynomial #in the form of a polynomial in (m,n,1/M,1/N) #For example, try: CanOpe(M*N-M-N,m,n,M,N); CanOpe:=proc(Ope,m,n,M,N) local gu,d1,d2,Ope1,i1,j1: Ope1:=numer(Ope): d1:=degree(Ope1,M): d2:=degree(Ope1,N): Ope1:=expand(subs({m=m-d1,n=n-d2},Ope1)/M^d1/N^d2): gu:=0: for i1 from ldegree(Ope1,M) to degree(Ope1,M) do for j1 from ldegree(Ope1,N) to degree(Ope1,N) do gu:=gu+factor(coeff(coeff(Ope1,M,i1),N,j1))*M^i1*N^j1: od: od: gu: end: #CanOpePos(Ope,m,n,M,N): Given a partial linear-recurrence operator with polynomial coefficients #divides it by an appropriate shiftp operaor to find an equivalent polynomial #in the form of a polynomial in (m,n,M,N) #For example, try: CanOpePos(M*N-M/m-N,m,n,M,N); CanOpePos:=proc(Ope,m,n,M,N) local gu,d1,d2,Ope1,i1,j1: Ope1:=numer(Ope): d1:=ldegree(Ope1,M): d2:=ldegree(Ope1,N): Ope1:=expand(subs({m=m-d1,n=n-d2},Ope1)/M^d1/N^d2): gu:=0: for i1 from 0 to degree(Ope1,M) do for j1 from 0 to degree(Ope1,N) do gu:=gu+factor(coeff(coeff(Ope1,M,i1),N,j1))*M^i1*N^j1: od: od: gu: end: #findrecR(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of rational numbers degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecR:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: 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: #FindrecR(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); FindrecR:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 1 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:=findrecR([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #Tp(f,x,p): The transformation T_p(f(x)) Tp:=proc(f,x,p) local i: normal( subs({seq(x[i]=p*x[i]/(1-(1-p)*x[i]),i=1..nops(x))},f)/ mul(1-(1-p)*x[i],i=1..nops(x))): end: #Sp(f,x,p): The transformation S_p(f(x)) Sp:=proc(f,x,p) local i: normal( subs({seq(x[i]=x[i]/(p+(1-p)*x[i]),i=1..nops(x))},f)/ mul(1+(1-p)*x[i]/p,i=1..nops(x))): end: #ApplyOper(Seq1,ope,n,N,i0): Given a sequence Seq1, applies the operator #ope(n,N) at Seq1[i] ApplyOper:=proc(Seq1,ope,n,N,i0) local i: add(Seq1[i0+i]*subs(n=i0,coeff(ope,N,i)),i=ldegree(ope,N)..degree(ope,N)): end: ez:=proc():print(`Mult1(ope1,ope2,n,N), Div1(ope1,ope2,n,N)`):end: #Div1(ope1,ope2,n,N): Given two operators ope1(n,N) and ope2(n,N) #outputs operators Q(n,N) and R(n,N) such that ope1=Q*ope2+R #and the order of R is less than the order of ope2 Div1:=proc(ope1,ope2,n,N) local d1,d2,i,ope1a,gu: d1:=degree(ope1,N): d2:=degree(ope2,N): if d1=n>0 (i.e. ballot walks). `): print(`By the way, the first few terms (starting at n=1), are`): print([op(1..10,gu)]): print(`b(n) satisfies the following linear recurrence `): print(add(coeff(ope,N,i)*b(n+i),i=ldegree(ope,N)..degree(ope,N))=0): print(`Skecth of proof: Analogous to the above.`): print(``): print(`This implies, via the Birkhoff-Trijinski-Proincare method, that the`): print(`the asymptotics of b(n) is: `): lu2:=AsyC(ope,n,N,k,[op(1..degree(ope,N),gu)],M): gu2a:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],M): print(lu2): print(`Disclaimer: The asymptotics is fully proved and rigorous, but`): print(`the constant in front is a non-rigorous estimate.`): print(`And in floating-point`, evalf(lu2)): print(`By dividing the asymptotic expression of b(n) by that of a(n)`): print(` we get the following`): print(`Corollary: The asymptotics of the probabibility of a walk `): print(`that winded up at (n,n) staying in m>=n is:`): LU:=normal(lu2/lu1): LU:=normal(subs(n=1/n,LU)): LU:=series(LU,n=0,k+1): LU:=add(coeff(LU,n,i)*n^i,i=0..k): LU:=subs(n=1/n,LU): print(LU+O(1/n^(k+1)) ): GU:=[seq(evalf(gu2a[i]/gu1a[i]/subs(n=i,LU)),i=nops(gu1a)-9..nops(gu1a))]: print(`----------------------------`): print(`Just as a check, the last 10 terms of the sequence of ratios of `): print(`the probabilities to the above asymptotics is before the`, M): print(`term is:`): print(GU): if abs(1-GU[nops(GU)])<1/100 then print(`as you can see, they are pretty close to 1`): else print(`this is rather disappointing, there must be something wrong.`): fi: print(`Disclaimer: The asymptotics is fully proved and rigorous, but`): print(`the constant in front is a non-rigorous estimate.`): print(`----------------------------`): print(`The whole thing took `, time()-t0, `seconds of CPU time.`): end: #LatticePathStory(Steps,Patience, s1): the story about Lattice Paths with Fundamental Steps #Steps: For example, try: LatticePathStory({[0,1],[1,0],[1,1]},10,15); LatticePathStory:=proc(Steps, Patience, s1) local gu,m,n,M,N,i,j,opeM,opeN,P,mu,ku,t0: t0:=time(): gu:=GH2((i,j)->LaW(Steps,[i,j]),m,n,M,N,Patience,s1) : if gu=FAIL or gu[1]=FAIL or gu[2]=FAIL then print(`Make the second and/or third arguments larger`): RETURN(FAIL): fi: opeM:=gu[1]: opeN:=gu[2]: print(`Theorem: Let F(m,n) be the number of lattice paths from (0,0) to (m,n)`): print(`in the positive quadrant of the two-dimensional square lattice`): print(`where the Fundamental Steps are`): print(Steps): print(`Let M and N be the shift operators in the m and n directions, respectively`): print(` (M f(m,n):=f(m+1,n), Nf(m,n):=f(m,n+1) ) `): print(` Then F(m,n) satisfies the following pure recurrences`): print(add(coeff(opeM,M,i)*F(m+i,n),i=ldegree(opeM,M)..degree(opeM,M))=0): print(``): print(add(coeff(opeN,N,i)*F(m,n+i),i=ldegree(opeN,N)..degree(opeN,N))=0): print(`Proof: We will only do the first recurrence`): print(`We have to prove that F(m,n) is annihilated by the operator, let's call it Q:=`): print(opeM): opeM:=numer(normal(opeM)): print(`or equivalenty`): opeM:=numer(normal(opeM)): print(opeM): print(`We know, by the obvious combinatorics, that F(m,n) is annihilated by`): P:=1-add(1/M^Steps[i][1]/N^Steps[i][2],i=1..nops(Steps)): P:=numer(normal(P)): print(P): mu:=CommSeq(P,opeM,m,n,M,N): print(`The sequence of successive commutators is`): print([op(2..nops(mu),mu)]): ku:=normal(mu[nops(mu)]/P): if degree(denom(ku),{M,N})<>0 then print(`Something is wrong with the proof`): else print(`As you can see, the last entry`, mu[nops(mu)], `is a multiple of `, P): print(`The proof follows by backwards induction, after checking the boundary conditions`): print(` QED .`): fi: print(`------------------------------------------------------`): print(`The whole thing took `, time()-t0, `seconds of CPU time`): end: #ReciPolStory(POL,x,y,Patience, s1): the story about the Maclaurin coefficients of the rational #function 1/POL, where P is a polynomial in x,y #For example, try: ReciPolStory(1-x-y,x,y,10,15); ReciPolStory:=proc(POL,x,y, Patience, s1) local gu,m,n,M,N,i,j,opeM,opeN,P,mu,ku,t0: t0:=time(): gu:=GH2Pol(POL,x,y,m,n,M,N,Patience,s1): if gu=FAIL or gu[1]=FAIL or gu[2]=FAIL then print(`Make the second and/or third arguments larger`): RETURN(FAIL): fi: opeM:=gu[1]: opeN:=gu[2]: print(`Theorem: Let F(m,n) be the coefficient of x^m*y^n in the Maclaurin series of`): print(1/POL): print(`F(m,n) satisfies the following pure recurrences`): print(add(coeff(opeM,M,i)*F(m+i,n),i=ldegree(opeM,M)..degree(opeM,M))=0): print(``): print(add(coeff(opeN,N,i)*F(m,n+i),i=ldegree(opeN,N)..degree(opeN,N))=0): print(`Proof: We will only do the first recurrence`): print(`Let M and N be the shift operators in the m and n directions, respectively`): print(` (M f(m,n):=f(m+1,n), Nf(m,n):=f(m,n+1) ) `): print(`We have to prove that F(m,n) is annihilated by the operator, let's call it Q:=`): print(opeM): opeM:=numer(normal(opeM)): print(`or equivalenty`): opeM:=numer(normal(opeM)): print(opeM): print(`We know, by the obvious algebra, that F(m,n) is annihilated by`): P:=subs({x=1/M,y=1/N},POL): P:=numer(normal(P)): print(P): mu:=CommSeq(P,opeM,m,n,M,N): print(`The sequence of successive commutators is`): print([op(2..nops(mu),mu)]): ku:=normal(mu[nops(mu)]/P): if degree(denom(ku),{M,N})<>0 then print(`Something is wrong with the proof`): else print(`As you can see, the last entry`, mu[nops(mu)], `is a multiple of `, P): print(`The proof follows by backwards induction, after checking the boundary conditions`): print(` QED .`): fi: print(`--------------------------------------------------------------------------------`): print(`The whole thing took `, time()-t0, `seconds of CPU time`): end: