###################################################################### ## QuarterPlane Save this file as QuarterPlane to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read QuarterPlane # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg at math dot rutgers dot edu. # ###################################################################### print(`First Written: April 17, 2007: tested for Maple 10 `): print(`Slightly revised June 22, 2021 `): print(): print(`This is QuarterPlane, A Maple package`): print(`To study walks on the Quarter Plane`): print(`accompanying Manuel Kauers and Doron Zeilberger's article: `): print(`"The Quasi-Holonomic Ansatz and Restricted Lattice Walks" `): print(): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/QuarterPlane .`): 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(`The supporting procedures are: Findrec, GenPureOper,`): print(` ApplyOper, GenOper, GuessOper1, GuessPureOper1 `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` QuarterPlane: A Maple package for studying `): print(`walks with arbitrary steps, confined to the quarter plane`): print(`The MAIN procedures are`): print(` OpeK, GuessRec0, W `): elif nargs=1 and args[1]=ApplyOper then print(`ApplyOper(n0,i0,j0,Steps,Ope,Sn,i,j): inputs n0,i0,j0, a set of steps Steps, and an operator Ope in Sn,i,j, applies it to the function W(n,i,j) at n0,i0,j0. Try:`): print(`ApplyOper(4,1,2,{[-1,0],[0,-1],[1,1]},OpeK(n,i,j,Sn),Sn,i,j);`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N,MaxC): `): print(`Given a list f tries to find a linear recurrence equation with`): print(`poly coffs. 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]=GenPureOper then print(`GenPureOper(n,a,b,N,d1): a recurrence operator `): print(`of the form P(n,a,b,N)`): print(`of degree d1`): print(`For example, try:GenPureOper(n,a,b,N,d1);`): elif nargs=1 and args[1]=GenOper then print(`GenOper(n,a,b,N,A,B,d1,d2): a generic recurrence operator `): print(`of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B)`): print(`with degree of P d1, and degrees of Q1 and Q2 d2`): print(`followed by the list of coefficients `): print(`For example, try:GenOper(n,a,b,N,A,B,1,1);`): elif nargs=1 and args[1]=GuessOper1 then print(`GuessOper1(n,a,b,N,A,B,d1,d2,Steps): `): print(`guesses a a recurrence operator `): print(`of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B)`): print(`with degree of P d1, and degrees of Q1 and Q2 d2`): print(`annihilating F(n,a,b):= The number of n-step walks `): print(`from the origin to [a,b] staying in`): print(`the quarter-plane using the steps in the set of Steps Steps`): print(`For example, try:`): print(`GuessOper1(n,a,b,N,A,B,1,1,{[0,1],[1,0],[0,-1],[-1,0]});`): elif nargs=1 and args[1]=GuessPureOper1 then print(`GuessPureOper1(n,a,b,N,d1,Steps): `): print(`guesses a a recurrence operator `): print(`of the form P(n,a,b,N)`): print(`of total degree d1 `): print(`annihilating F(n,a,b):= The number of n-step walks `): print(`from the origin to [a,b] staying in`): print(`the quarter-plane using the steps in the set of Steps Steps`): print(`For example, try:`): print(`GuessPureOper1(n,a,b,N,1,{[0,1],[1,0],[0,-1],[-1,0]});`): elif nargs=1 and args[1]=GuessRec0 then print(`GuessRec0(Steps,a,b,n,N,MaxC): guesses a recurrence operator`): print(`ope(n,N) of complexity <=MaxC`): print(`for the sequence: number of ways of going from the origin`): print(`to [a,b] using data for n<=L, for example, try:`): print(`GuessRec0({[1,1],[-1,0],[0,-1]},0,0,n,N,4);`): elif nargs=1 and args[1]=OpeK then print(`OpeK(n,i,j,Sn): the precomuted pure operator for Kreweras walks, Type:`): print(`OpeK(n,i,j,Sn);`): elif nargs=1 and args[1]=W then print(`W(n,a,b,Steps): the number of n-step walks in the positive `): print(`discrete (square-lattice) quarter-plane staring at the`): print(`origin and ending at [a,b] using the set of steps in Steps`): print(`For example, try: W(3,2,1,{[1,0],[0,1]});`): print(`For Kreweras, try: W(9,0,0,{[-1,0],[0,-1],[1,1]});`): print(`For Gessel, try: W(8,0,0,{[-1,0],[1,0],[-1,-1],[1,1]});`): else print(`There is no such thing as`, args): fi: end: #Start pre-computed #OpeK(n,i,j,Sn): the precomuted pure operator for Kreweras walks OpeK:=proc(n,i,j,Sn): 4251528*(1+n)*(2+n)*(3+n)*(4+n)*(5+n)*(6+n)*(7+n)*(8+n)*(9+n)*(8*i^3-12*i^2*j+10*i^2*n-12*i*j^2-10*i*j*n+8*j^3+10*j^2*n+100*i^2-136*i*j+10*i*n+100*j^2+10*j*n+35*n^2+76*i+76*j+785*n+4384)+78732*(4+n)*(5+n)*(6+n)*(7+n)*(8+n)*(9+n)*(96*i^6-288*i^5*j+264*i^5*n-72*i^4*j^2-660*i ^4*j*n+180*i^4*n^2+624*i^3*j^3+264*i^3*j^2*n-360*i^3*j*n^2-64*i^3*n^3-72*i^2*j^4+264*i^2*j^3*n+540*i^2*j^2*n^2+96*i^2*j*n^3-80*i^2*n^4-288*i*j^5-660*i*j^4*n-360*i*j^3*n^2+96*i*j^2*n^3+80*i*j*n^4+96*j^6+264*j^5*n+180*j^4*n^2-64*j^3*n^3-80*j^2*n^4+2088*i^5-6084*i^4*j+3144*i^ 4*n+3384*i^3*j^2-7080*i^3*j*n-540*i^3*n^2+3384*i^2*j^3+9036*i^2*j^2*n+1350*i^2*j*n^2-1760*i^2*n^3-6084*i*j^4-7080*i*j^3*n+1350*i*j^2*n^2+2048*i*j*n^3-80*i*n^4+2088*j^5+3144*j^4*n-540*j^3*n^2-1760*j^2*n^3-80*j*n^4-280*n^5+12312*i^4-31752*i^3*j+3868*i^3*n+39204*i^2*j^2+462*i ^2*j*n-9700*i^2*n^2-31752*i*j^3+462*i*j^2*n+14290*i*j*n^2-1568*i*n^3+12312*j^4+3868*j^3*n-9700*j^2*n^2-1568*j*n^3-11635*n^4+36624*i^3-38736*i^2*j+14750*i^2*n-38736*i*j^2-3536*i*j*n-7180*i*n^2+36624*j^3+14750*j^2*n-7180*j*n^2-192787*n^3+157674*i^2-237270*i*j+15962*i*n+ 157674*j^2+15962*j*n-1609205*n^2+88410*i+88410*j-6839845*n-11889960)*Sn^3+4374*(7+n)*(8+n)*(9+n)*(128*i^9-576*i^8*j+544*i^8*n+288*i^7*j^2-2176*i^7*j*n+768*i^7*n^2+1680*i^6*j^3+1768*i^6*j^2*n-2688*i^6*j*n^2+168*i^6*n^3-1584*i^5*j^4+2312*i^5*j^3*n+3456*i^5*j^2*n^2-504*i^5*j* n^3-528*i^5*n^4-1584*i^4*j^5-4352*i^4*j^4*n-1920*i^4*j^3*n^2+2304*i^4*j^2*n^3+1320*i^4*j*n^4-360*i^4*n^5+1680*i^3*j^6+2312*i^3*j^5*n-1920*i^3*j^4*n^2-3768*i^3*j^3*n^3-528*i^3*j^2*n^4+720*i^3*j*n^5+64*i^3*n^6+288*i^2*j^7+1768*i^2*j^6*n+3456*i^2*j^5*n^2+2304*i^2*j^4*n^3-528* i^2*j^3*n^4-1080*i^2*j^2*n^5-96*i^2*j*n^6+80*i^2*n^7-576*i*j^8-2176*i*j^7*n-2688*i*j^6*n^2-504*i*j^5*n^3+1320*i*j^4*n^4+720*i*j^3*n^5-96*i*j^2*n^6-80*i*j*n^7+128*j^9+544*j^8*n+768*j^7*n^2+168*j^6*n^3-528*j^5*n^4-360*j^4*n^5+64*j^3*n^6+80*j^2*n^7+4544*i^8-19904*i^7*j+13360* i^7*n+19952*i^6*j^2-50840*i^6*j*n+6296*i^6*n^2+20608*i^5*j^3+67872*i^5*j^2*n-20040*i^5*j*n^2-14984*i^5*n^3-47584*i^4*j^4-33808*i^4*j^3*n+61788*i^4*j^2*n^2+40808*i^4*j*n^3-14448*i^4*n^4+20608*i^3*j^5-33808*i^3*j^4*n-97856*i^3*j^3*n^2-17576*i^3*j^2*n^3+30480*i^3*j*n^4+1656*i ^3*n^5+19952*i^2*j^6+67872*i^2*j^5*n+61788*i^2*j^4*n^2-17576*i^2*j^3*n^3-42552*i^2*j^2*n^4-3564*i^2*j*n^5+3440*i^2*n^6-19904*i*j^7-50840*i*j^6*n-20040*i*j^5*n^2+40808*i*j^4*n^3+30480*i*j^3*n^4-3564*i*j^2*n^5-3728*i*j*n^6+80*i*n^7+4544*j^8+13360*j^7*n+6296*j^6*n^2-14984*j^5 *n^3-14448*j^4*n^4+1656*j^3*n^5+3440*j^2*n^6+80*j*n^7+280*n^8+56240*i^7-232648*i^6*j+66664*i^6*n+335112*i^5*j^2-229008*i^5*j*n-153952*i^5*n^2-165176*i^4*j^3+561804*i^4*j^2*n+455476*i^4*j*n^2-225902*i^4*n^3-165176*i^3*j^4-836224*i^3*j^3*n-214648*i^3*j^2*n^2+504964*i^3*j*n^3 +4952*i^3*n^4+335112*i^2*j^5+561804*i^2*j^4*n-214648*i^2*j^3*n^2-667650*i^2*j^2*n^3-44436*i^2*j*n^4+52160*i^2*n^5-232648*i*j^6-229008*i*j^5*n+455476*i*j^4*n^2+504964*i*j^3*n^3-44436*i*j^2*n^4-63932*i*j*n^5+3248*i*n^6+56240*j^7+66664*j^6*n-153952*j^5*n^2-225902*j^4*n^3+4952 *j^3*n^4+52160*j^2*n^5+3248*j*n^6+19510*n^7+218000*i^6-804696*i^5*j-671048*i^5*n+1713288*i^4*j^2+2148824*i^4*j*n-1721950*i^4*n^2-2366888*i^3*j^3-1104560*i^3*j^2*n+4084028*i^3*j*n^2-244996*i^3*n^3+1713288*i^2*j^4-1104560*i^2*j^3*n-5209314*i^2*j^2*n^2-127704*i^2*j*n^3+209995 *i^2*n^4-804696*i*j^5+2148824*i*j^4*n+4084028*i*j^3*n^2-127704*i*j^2*n^3-387439*i*j*n^4+45392*i*n^5+218000*j^6-671048*j^5*n-1721950*j^4*n^2-244996*j^3*n^3+209995*j^2*n^4+45392*j*n^5+596967*n^6-1030256*i^5+3521288*i^4*j-6409640*i^4*n-1913480*i^3*j^2+16087168*i^3*j*n-3255704 *i^3*n^2-1913480*i^2*j^3-20141940*i^2*j^2*n+1642854*i^2*j*n^2-2944298*i^2*n^3+3521288*i*j^4+16087168*i*j^3*n+1642854*i*j^2*n^2+1859816*i*j*n^3+128707*i*n^4-1030256*j^5-6409640*j^4*n-3255704*j^3*n^2-2944298*j^2*n^3+128707*j*n^4+10489063*n^5-9378548*i^4+24654784*i^3*j-\ 16287660*i^3*n-30540900*i^2*j^2+13967190*i^2*j*n-40622871*i^2*n^2+24654784*i*j^3+13967190*i*j^2*n+40141659*i*j*n^2-3170468*i*n^3-9378548*j^4-16287660*j^3*n-40622871*j^2*n^2-3170468*j*n^3+115857999*n^4-30133640*i^3+31732608*i^2*j-193291854*i^2*n+31732608*i*j^2+214924128*i*j *n-37672845*i*n^2-30133640*j^3-193291854*j^2*n-37672845*j*n^2+824266699*n^3-338836500*i^2+404147016*i*j-168223842*i*n-338836500*j^2-168223842*j*n+3689510562*n^2-281471472*i-281471472*j+9499266216*n+10767717504)*Sn^6+27*(-9+2*i-j-n)*(9+i-2*j+n)*(12+i+j+n)*(256*i^9-1152*i^8* j+1088*i^8*n+576*i^7*j^2-4352*i^7*j*n+1536*i^7*n^2+3360*i^6*j^3+3536*i^6*j^2*n-5376*i^6*j*n^2+80*i^6*n^3-3168*i^5*j^4+4624*i^5*j^3*n+6912*i^5*j^2*n^2-240*i^5*j*n^3-1760*i^5*n^4-3168*i^4*j^5-8704*i^4*j^4*n-3840*i^4*j^3*n^2+4800*i^4*j^2*n^3+4400*i^4*j*n^4-1200*i^4*n^5+3360*i ^3*j^6+4624*i^3*j^5*n-3840*i^3*j^4*n^2-9200*i^3*j^3*n^3-1760*i^3*j^2*n^4+2400*i^3*j*n^5+256*i^3*n^6+576*i^2*j^7+3536*i^2*j^6*n+6912*i^2*j^5*n^2+4800*i^2*j^4*n^3-1760*i^2*j^3*n^4-3600*i^2*j^2*n^5-384*i^2*j*n^6+320*i^2*n^7-1152*i*j^8-4352*i*j^7*n-5376*i*j^6*n^2-240*i*j^5*n^3 +4400*i*j^4*n^4+2400*i*j^3*n^5-384*i*j^2*n^6-320*i*j*n^7+256*j^9+1088*j^8*n+1536*j^7*n^2+80*j^6*n^3-1760*j^5*n^4-1200*j^4*n^5+256*j^3*n^6+320*j^2*n^7+9536*i^8-41600*i^7*j+27872*i^7*n+41360*i^6*j^2-105712*i^6*j*n+4288*i^6*n^2+43120*i^5*j^3+140928*i^5*j^2*n-15168*i^5*j*n^2-\ 61760*i^5*n^3-98752*i^4*j^4-70496*i^4*j^3*n+134664*i^4*j^2*n^2+163400*i^4*j*n^3-55760*i^4*n^4+43120*i^3*j^5-70496*i^3*j^4*n-259408*i^3*j^3*n^2-70400*i^3*j^2*n^3+116800*i^3*j*n^4+9872*i^3*n^5+41360*i^2*j^6+140928*i^2*j^5*n+134664*i^2*j^4*n^2-70400*i^2*j^3*n^3-164640*i^2*j^2 *n^4-18408*i^2*j*n^5+17120*i^2*n^6-41600*i*j^7-105712*i*j^6*n-15168*i*j^5*n^2+163400*i*j^4*n^3+116800*i*j^3*n^4-18408*i*j^2*n^5-18272*i*j*n^6+320*i*n^7+9536*j^8+27872*j^7*n+4288*j^6*n^2-61760*j^5*n^3-55760*j^4*n^4+9872*j^3*n^5+17120*j^2*n^6+320*j*n^7+1120*n^8+122912*i^7-\ 505168*i^6*j+47824*i^6*n+723552*i^5*j^2-203232*i^5*j*n-821616*i^5*n^2-356768*i^4*j^3+1280616*i^4*j^2*n+2279688*i^4*j*n^2-1031060*i^4*n^3-356768*i^3*j^4-2422624*i^3*j^3*n-1064832*i^3*j^2*n^2+2266120*i^3*j*n^3+138920*i^3*n^4+723552*i^2*j^5+1280616*i^2*j^4*n-1064832*i^2*j^3*n ^2-3039780*i^2*j^2*n^3-354540*i^2*j*n^4+369660*i^2*n^5-505168*i*j^6-203232*i*j^5*n+2279688*i*j^4*n^2+2266120*i*j^3*n^3-354540*i*j^2*n^4-428484*i*j*n^5+16352*i*n^6+122912*j^7+47824*j^6*n-821616*j^5*n^2-1031060*j^4*n^3+138920*j^3*n^4+369660*j^2*n^5+16352*j*n^6+86020*n^7+ 143264*i^6-754224*i^5*j-4932720*i^5*n+4090512*i^4*j^2+14220864*i^4*j*n-9502848*i^4*n^2-7558880*i^3*j^3-7158816*i^3*j^2*n+21921264*i^3*j*n^2+764920*i^3*n^3+4090512*i^2*j^4-7158816*i^2*j^3*n-28335744*i^2*j^2*n^2-3505560*i^2*j*n^3+4000500*i^2*n^4-754224*i*j^5+14220864*i*j^4*n +21921264*i*j^3*n^2-3505560*i*j^2*n^3-5234040*i*j*n^4+334044*i*n^5+143264*j^6-4932720*j^5*n-9502848*j^4*n^2+764920*j^3*n^3+4000500*j^2*n^4+334044*j*n^5+2897468*n^6-11321472*i^5+33597120*i^4*j-43718268*i^4*n-18009696*i^3*j^2+105838584*i^3*j*n-34992*i^3*n^2-18009696*i^2*j^3-\ 133285356*i^2*j^2*n-18824616*i^2*j*n^2+20703300*i^2*n^3+33597120*i*j^4+105838584*i*j^3*n-18824616*i*j^2*n^2-34398960*i*j*n^3+3324300*i*n^4-11321472*j^5-43718268*j^4*n-34992*j^3*n^2+20703300*j^2*n^3+3324300*j*n^4+55920771*n^5-80397744*i^4+204321072*i^3*j-14473344*i^3*n-\ 252690336*i^2*j^2-52970364*i^2*j*n+21490272*i^2*n^2+204321072*i*j^3-52970364*i*j^2*n-107422416*i*j*n^2+13934280*i*n^3-80397744*j^4-14473344*j^3*n+21490272*j^2*n^2+13934280*j*n^3+676455675*n^4-34297872*i^3-64163736*i^2*j-230206104*i^2*n-64163736*i*j^2-63941616*i*j*n-\ 16857216*i*n^2-34297872*j^3-230206104*j^2*n-16857216*j*n^2+5251873575*n^3-706741884*i^2+267970932*i*j-351219204*i*n-706741884*j^2-351219204*j*n+25552376145*n^2-882810684*i-882810684*j+71210374074*n+86992033968)*Sn^9-(-21+4*i-2*j-2*n)*(-12+2*i-j-n)*(-9+2*i-j-n)*(9+i-2*j+n)* (12+i-2*j+n)*(12+i+j+n)*(15+i+j+n)*(21+2*i-4*j+2*n)*(27+2*i+2*j+2*n)*(8*i^3-12*i^2*j+10*i^2*n-12*i*j^2-10*i*j*n+8*j^3+10*j^2*n+70*i^2-106*i*j+10*i*n+70*j^2+10*j*n+35*n^2+46*i+46*j+575*n+2344)*Sn^12: end: #End pre-computed ################Start Findrec with(linalg): #with(SolveTools); #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: #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,i,n1: 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 if {seq(expand(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..ORDER)), n1=1..nops(f)-ORDER)}={0} then RETURN(ope): else print({seq(expand(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..ORDER)), n1=1..nops(f)-ORDER)}): fi: fi: od: od: FAIL: end: ################End Findrec and its offsprings #W(n,a,b,Steps): the number of n-step walks in the positive #discrete (square-lattice) quarter-plane staring at the #origin and ending at [a,b] using the set of steps in Steps #For example, try: W(3,2,1,{[1,0],[0,1]}); W:=proc(n,a,b,Steps) local i: option remember: if n<0 then RETURN(0): fi: if n=0 then if a=0 and b=0 then RETURN(1): else RETURN(0): fi: fi: if a<0 or b<0 then RETURN(0): fi: add(W(n-1,a-Steps[i][1],b-Steps[i][2],Steps),i=1..nops(Steps)): end: #GenOper(n,a,b,N,A,B,d1,d2): a recurrence operator #of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B) #with degree of P d1, and degrees of Q1 and Q2 d2 #For example, try:GenOper(n,a,b,N,A,B,1,1); GenOper:=proc(n,a,b,N,A,B,d1,d2) local ope,var,i1,i2,i3,i4,i5,i6,c1,c2,c3: ope:=0: var:={}: for i1 from 0 to d1 do for i2 from 0 to d1-i1 do ope:=ope+c1[i1,i2]*n^i1*N^i2: var:=var union {c1[i1,i2]}: od: od: for i1 from 0 to d2 do for i2 from 0 to d2-i1 do for i3 from 0 to d2-i1-i2 do for i4 from 0 to d2-i1-i2-i3 do for i5 from 0 to d2-i1-i2-i3-i4 do for i6 from 0 to d2-i1-i2-i3-i4-i5 do ope:=ope+a*c2[i1,i2,i3,i4,i5,i6]*n^i1*N^i2*a^i3*b^i4*A^i5*B^i6: var:=var union {c2[i1,i2,i3,i4,i5,i6]}: od: od: od: od: od: od: for i1 from 0 to d2 do for i2 from 0 to d2-i1 do for i3 from 0 to d2-i1-i2 do for i4 from 0 to d2-i1-i2-i3 do for i5 from 0 to d2-i1-i2-i3-i4 do for i6 from 0 to d2-i1-i2-i3-i4-i5 do ope:=ope+b*c3[i1,i2,i3,i4,i5,i6]*n^i1*N^i2*a^i3*b^i4*A^i5*B^i6: var:=var union {c3[i1,i2,i3,i4,i5,i6]}: od: od: od: od: od: od: ope,var: end: #GuessOper1(n,a,b,N,A,B,d1,d2,Steps): #guesses a a recurrence operator #of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B) #with degree of P d1, and degrees of Q1 and Q2 d2 #annihilating F(n,a,b):= The number of n-step walks #from the origin to [a,b] staying in #the quarter-plane using the steps in the set of Steps Steps #For example, try:GuessOper1(n,a,b,N,A,B,1,1,{[0,1],[1,0],[0,-1],[-1,0]}); GuessOper1:=proc(n,a,b,N,A,B,d1,d2,Steps) local ope,var,n1,a1,b1,eq,eq1,K,j,mu,dega,degb,degn,coe,var1: ope:=GenOper(n,a,b,N,A,B,d1,d2): var:=ope[2]: ope:=ope[1]: eq:={}: for K from 0 while nops(eq)0 then gu:=gu union {gu1[i]}: fi: od: gu: end: #GuessRec0(Steps,a,b,n,N,MaxC): guesses a recurrence operator #ope(n,N), of complexity MaxC #for the sequence: number of ways of going from the origin #to [a,b] using data for n<=L, for example, try: #GuessRec0({[1,1],[-1,0],[0,-1]},0,0,n,N,40); GuessRec0:=proc(Steps,a,b,n,N,MaxC) local gu,i,n1,a1,b1,r,j,L,ope: L:=(2+MaxC)^2+4: gu:=[seq(W(n1,a,b,Steps),n1=0..L)]: for i from 1 while gu[i]=0 do od: a1:=i: gu:=[op(a1..nops(gu),gu)]: for i from 2 to nops(gu) while gu[i]=0 do od: b1:=i-1: if {seq(seq(gu[j+b1*r],r=0..trunc((nops(gu)-j)/b1)), j=2..b1)}<>{0} then RETURN(FAIL): fi: L:=a1+b1*L: gu:=[seq(W(n1,a,b,Steps),n1=0..L)]: gu:=[seq(gu[a1+b1*r],r=0..trunc((nops(gu)-a1-1)/b1))]: ope:=Findrec(gu,n,N,MaxC): if ope=FAIL then RETURN(FAIL): else RETURN(ope,a1-1,b1): fi: end: #GenSP2(x,y,d,a): A generic symmettic polynomial in (x,y) of total degree d with coefficients parametrized by a. It also outouts the set of coefficients GenSP2:=proc(x,y,d,a) end: #ApplyOper(n0,i0,j0,Steps,Ope,Sn,i,j): inputs n0,i0,j0, a set of steps Steps, and an operator Ope in Sn,i,j, applies it to the function W(n,i,j) at n0,i0,j0. Try: #ApplyOper(4,1,2,{[-1,0],[0,-1],[1,1]},OpeK(n,i,j,Sn),Sn,i,j); ApplyOper:=proc(n0,i0,j0,Steps,Ope,Sn,i,j) local a1: add(subs({i=i0,j=j0,n=n0},coeff(Ope,Sn,a1))*W(n0+a1,i0,j0,Steps),a1=0..degree(Ope,Sn)): end: