##################################################################### ## W2D Save this file as W2D to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read W2D # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg at math dot rutgers dot edu. # ###################################################################### with(combinat): print(`First Written: Jan. , 2015: tested for Maple 16 `): print(`Version of Jan. 13, 2015: `): print(): print(`This is W2D, A Maple package`): print(`for the enumeration and probability of 2 dimensional walks restricted to the`): print(`quarter plane`): print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's article: `): print(`"The Method(!) of Guess and Check" `): print(): print(`The most current version is available at:`): print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/W2D .`): 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(): Digits:=40: ezra1:=proc() if args=NULL then print(` The supporting procedures are:`): print(` EmpirF, NuWalks2D, OneStepFeK2D, PrintFE2D, Seq2DxySlow, TikvaLight`): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` W2D: A Maple package for the enumeration and probability of`): print(`two-dimensional walks restricted to the positive quarter-plane `): print(`The MAIN procedures are: FirstK2D, InfoA2D, InfoH2D, InfoH2Dv, InfoH2DvWP, Sefer, Seq2D, Seq2Dxy, Tikva, Wfe2D `): print(``): elif nargs=1 and args[1]=FirstK2D then print(` FirstK2D(F,t,x,y,ftxy,C,Cx,Cy,H,K): Given `): print(` a functional equation F phrased in terms of C[i,j], i=0..H[1], j=0.., H[2], and Cx[i], i=0..H[1], and Cy[i], i=0..H[2] `): print(` finds the first K terms of the formal power series in t, with polynomials in x and y, that satisfies it. Try: `): print(` lu:=Wfe2D({[1,1],[-1,0],[0,-1]},ftxy, t,x,y,Cx,Cy,C); `): print(` FirstK2D(lu[1],t,x,y,ftxy,C,Cx,Cy,lu[2],10); `): print(`Also try: `): print(` lu:=Wfe2D({[1,1],[-1,-1],[1,0],[-1,0]},ftxy, t,x,y,Cx,Cy,C); `): print(` FirstK2D(lu[1],t,x,y,ftxy,C,Cx,Cy,lu[2],10); `): elif nargs=1 and args[1]=InfoA2D then print(` InfoA2D(S,t,P,K): `): print(`inputs a set of steps, symbols t and P (t is the variable, and P is the generating function) `): print(` and integer K, outputs, if possible, the pair [P0,P1], where P0 is the guessed algebraic equation `): print(`obeyed by the generating function of the the sequence enumerating walks that start at [0,0], use the `): print(` step-set S, and return to [0,0], all the while staying in the first quadrant, and `): print(` P1 is the algebraic equation of the generating function of the enumerating`): print(`sequence for such walks that end anywhere you want (in the first quardrant, of course). `): print(`If it returns FAIL it does not (necessarily) mean that there are no such equations, only that K terms do not suffice`): print(`to guess it. `): print(` For example, For the Kreweras walks try: `): print(` InfoA2D({[1,1],[-1,0],[0,-1]},t,P,60); `): print(` For reverse Kreweras walks try: `): print(` InfoA2D({[-1,-1],[1,0],[0,1]},t,P,60); `): print(` For the Gessel walks, try: `): print(`InfoA2D({[1,1],[-1,-1],[1,0],[-1,0]},t,P,60); `): elif nargs=1 and args[1]=InfoH2D then print(` InfoH2D(S,n,N,MaxC): `): print(`inputs a set of steps, symbols n and N (N is the shift operator in n, i.e. Na(n):=a(n+1)), and a positive `): print(` integer MaxC, outputs, if possible, the pair [Ope0,Ope1], where Ope0 is the holonomic representation `): print(`with ORDER+Degree<=MaxC of the sequence enumerating walks that start at [0,0], use the `): print(` step-set S, and return to [0,0], all the while staying in the first quadrant, and `): print(` Ope1 is the holonomic representation of such walks that end anywhere you want (in the first quardrant, of course). `): print(` For example, For the Kreweras walks try: `): print(` InfoH2D({[1,1],[-1,0],[0,-1]},n,N,10); `): print(` For reverse Kreweras walks try: `): print(` InfoH2D({[-1,-1],[1,0],[0,1]},n,N,11); `): print(` For the Gessel walks, try: `): print(`InfoH2D({[1,1],[-1,-1],[1,0],[-1,0]},n,N,5); `): elif nargs=1 and args[1]=InfoH2Dv then print(`InfoH2Dv(S,n,MaxC): Verbose version of InfoH2D(S,n,N,MaxC) (q.v.) `): print(` For example, For the Kreweras walks try: `): print(` InfoH2Dv({[1,1],[-1,0],[0,-1]},n,10); `): print(` For reverse Kreweras walks try: `): print(` InfoH2Dv({[-1,-1],[1,0],[0,1]},n,11); `): print(` For the Gessel walks, try: `): print(` InfoH2Dv({[1,1],[-1,-1],[1,0],[-1,0]},n,5); `): elif nargs=1 and args[1]=InfoH2DvWP then print(`InfoH2DvWP(S,n,MaxC): Verbose version of InfoH2D(S,n,N,MaxC) (q.v.) with a sketch of proof`): print(`For example, For the Kreweras walks try:`): print(`InfoH2DvWP({[1,1],[-1,0],[0,-1]},n,10); `): elif nargs=1 and args[1]=NuWalks2D then print(`NuWalks2D(Steps,t,x): Given a set of 2D steps Steps (pairs of positive or negative `): print(`integers or 0), a non-neg. integer t and a pair x returns the number of walks`): print(`(Sequences in Steps) after time t, arrive at x, and always stay in the non-negative quarter-plane`): print(`For example, try:`): print(` NuWalks2D({[1,1],[-1,0],[0,-1]},12,[0,0]); `): elif nargs=1 and args[1]=OneStepFeK2D then print(`OneStepFeK2D(g,F,t,x,y,ftxy,C,Cx,Cy,H,K): given a formal power series g in t, x, and y `): print(` and a functional equation F phrased in terms of C[i], i=0..H, `): print(`and a positive integer K, applies one step of the equation`): print(` g->F(g) and truncates it to K terms in t. Try: `): print(` lu:=Wfe2D({[1,0],[-1,0],[0,1],[0,-1]},ftxy, t,x,y,Cx,Cy,C);`): print(` OneStepFeK2D(1,lu[1],t,x,y,ftxy,C,Cx,Cy,lu[2],4); `): elif nargs=1 and args[1]=PrintFE2D then print(` PrintFE2D(eq,ftxy,t,x,y,C,Cx,Cy,P):prints nicely the functional equation eq phrased in terms of ftx,t,x,C. `): print(` Try: `): print(` PrintFE2D(Wfe2D({[1,0],[0,1],[-1,-1]},ftxy, t,x,y,Cx,Cy,C),ftxy,t,x,y,C,Cx,Cy,P,Coeff); `): elif nargs=1 and args[1]=PrintFE2Ds then print(` PrintFE2Ds(eq,ftxy,t,x,y,C,Cx,Cy,P):prints nicely the functional equation eq phrased in terms of ftx,t,x,C. `): print(`For the case where the set of steps S consists of short steps. `): print(` Try: `): print(` PrintFE2Ds(Wfe2D({[1,0],[0,1],[-1,-1]},ftxy, t,x,y,Cx,Cy,C),ftxy,t,x,y,C,Cx,Cy,P,Coeff); `): elif nargs=1 and args[1]=Sefer then print(`Sefer(k,MaxC): inputs a positive integer k between 3 and 8 and outputs all the theorems`): print(`about walks with k small steps, where the steps are subsets of `): print(` {[1,0],[-1,0],[0,1],[0,-1],[1,1],[1,-1],[-1,1],[-1,-1]}. Try `): print(` Sefer(3,10);`): elif nargs=1 and args[1]=Seq2D then print(` Seq2D(Steps,N,e): the enumerating sequence for the number of walks with i steps drawn from the set of steps Steps, `): print(` starting at [0,0] and ending at the point e, and always staying in the non-negative quarter-plane `): print(` Seq2D({[1,1],[-1,-1],[1,0],[-1,0]},20,[0,0]); `): print(` Seq2D({[1,1],[-1,0],[0,-1]},21,[0,0]); `): elif nargs=1 and args[1]=Seq2Dxy then print(`Seq2Dxy(S,N,x,y): inputs a set of steps and finds the first N terms of the enumerating polynomials for`): print(`n-step walks in the non-negative quarter plane, according to the weight x^dest[1]*y^dest[2], where`): print(`dest is the destination`): print(`Try:`): print(`Seq2Dxy({[1,1],[-1,-1],[1,0],[-1,0]},20,x,y);`): print(`Seq2Dxy({[1,1],[-1,0],[0,-1]},21,x,y);`): elif nargs=1 and args[1]=Seq2DxySlow then print(`Seq2DxySlow(S,N,x,y): A slow version of Seq2DxySlow (q.v.), using the Functional Equation.`): print(`Try:`): print(`Seq2DxySlow({[1,1],[-1,-1],[1,0],[-1,0]},20,x,y);`): print(`Seq2DxySlow({[1,1],[-1,0],[0,-1]},21,x,y);`): elif nargs=1 and args[1]=Tikva then print(`Tikva(S,MaxC): inputs a set of 2-dimensional steps, S, a positive integer MaxC`): print(`finds out whether the generating function of the first-quadrant walks f_S(x,y) is such that`): print(`f_S(x,0) and f_S(0,y) are holonomic in x,y. It returns the [[DegreeX,OrderX],[DegreY,OrderY]] `): print(`or FAIL `); print(`Try:`): print(`Tikva({[1,0],[0,1],[-1,0],[0,-1]},9);`): print(`Also try:`): print(`Tikva({[-1,0],[0,-1],[1,1]},10);`): elif nargs=1 and args[1]=TikvaLight then print(`TikvaLight(S,MaxC): inputs a set of 2-dimensional steps, S, a positive integer MaxC`): print(`finds out whether the generating function of the first-quadrant walks f_S(x,y) is such that`): print(`f_S(x,0) is holonomic in x Not as colclusive as Tivka(S,MaxC). It returns [Degree, Order]`): print(`or FAIL `); print(`Try:`): print(`TikvaLight({[1,0],[0,1],[-1,0],[0,-1]},9);`): print(`Also try:`): print(`TikvaLight({[-1,0],[0,-1],[1,1]},10);`): elif nargs=1 and args[1]=Wfe2D then print(` Wfe2D(S,ftxy, t,x,y,Cx,Cy,C): inputs a set of pairs of integers S, `): print(` indicating the allowed steps, `): print(` a symbol ftxy (denoting f(t,x,y), see below), variables t, and x, and y `): print(` symbols Cx,Cy, C `): print(` where Cx[i] denotes the coefficient of x^i in f(t,x,y) and `): print(`Cy[j] denotes the coefficient of y^i in f(t,x,y) and `): print(` and `): print(` C[i,j] denotes the coefficient of x^i*y^j `): print(` outputs the pair [eq,[H1,H2]], `): print(` where eq is functional equation for f(t,x,y), the generating `): print(` function of walks where the coefficient of t^i *x^j1*y^j2 is the number of ways `): print(` of walking i steps on the positive quarter-plane starting at the origin `): print(` and winding up at [j1,j2] and never going to outside the non-negative quater-plane `): print(` C[j1,j2] denotes the coefficient of x^j1*y^j2 in f(t,x,y), `): print(` H1 is the highest j1 that shows up as a C[j1,-] and H2 is the highest j2 that shows up a a C[-,j2] `): print(` Try: `): print(` Wfe2D({[1,0],[-1,0],[0,1],[0,-1]},ftxy, t,x,y,Cx,Cy,C); `): print(` Wfe2D({[1,1],[-1,-1],[1,0],[-1,0]},ftxy, t,x,y,Cx,Cy,C); `): print(` Wfe2D({[1,1],[-1,0],[1,0]},ftxy, t,x,y,Cx,Cy,C); `): else print(`There is no such thing as`, args): fi: end: sn:=proc(resh,n1): -1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2): end: #Zinn(resh): Zinn-Justin's method to estimate #the C,mu, and theta such that #resh[i] is appx. Const*mu^i*i^theta #For example, try: #Zinn([seq(5*i*2^i,i=1..30)]); Zinn:=proc(resh) local s1,s2,theta,mu,n1,i: if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then RETURN(FAIL): fi: n1:=nops(resh)-1: s1:=sn(resh,n1): s2:=sn(resh,n1-1): theta:=evalf(2*(s1+s2)/(s1-s2)^2): mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): [theta,mu]: end: ###############From Findrec with(linalg): #SeqFromRec(ope,n,N,Ini,K): given ope=[operator,Ini] #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,K) local ope1,gu,L,n1,j1,Ini: Ini:=ope[2]: ope1:=Yafe(ope[1],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: #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 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 if {seq(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..degree(ope,N)),n1=7..nops(f)-degree(ope,N))}={0} then RETURN([ope, [op(1..degree(ope,N),f)]]): fi: fi: od: od: FAIL: end: #KamaArokh(C): the length of the sequence needed to guess a recurrence of complexity C KamaArokh:=proc(C) local DEGREE: max(seq((2+DEGREE)*(1+C-DEGREE)+5,DEGREE=0..C)): end: ###############End from Findrec ####guess polynomial Di:=proc(L) local i: [ seq( L[i+1]-L[i],i=1..nops(L)-1)]: end: #GP(L,x): The unique polynomial P(x) #such that P(i)=L[i+1] for i from 0 to nops(L)-1 GP:=proc(L,x) local i,L1,P: L1:=L: P:=0: for i from 1 to nops(L) do P:=P+L1[1]*binomial(x,i-1): L1:=Di(L1): if convert(L1,set)={0} then RETURN(factor(expand(P))): fi: od: FAIL: end: ###End from GP #empir(gu,degx,degP,x,P) #to "fit" an algebraic equation F(P(x),x)=0 of degree #degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i empir:=proc(gu,degx,degP,x,P) local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip,var2,mu1,halev,vu,i,vu1: if (1+degx)*(1+degP) > nops(gu)-3 then RETURN(`sequence too small`): fi: F:=0: var:={}: for i1 from 0 to degx do for i2 from 0 to degP do F:=F+a[i1,i2]*x^i1*P^i2: var:=var union {a[i1,i2]}: od: od: cand:=0: for i1 from 0 to nops(gu)-1 do cand:=cand+op(i1+1,gu)*x^i1: od: lu:=subs(P=cand,F): lu:=taylor(lu,x=0,nops(gu)-1): eq:={}: for i1 from 0 to nops(gu)-2 do eq:=eq union {coeff(lu,x,i1)=0} od: mu:=solve(eq,var): F:=subs(mu,F): if F=0 then RETURN(FAIL): fi: var2:={}: for mu1 in mu do if op(1,mu1)=op(2,mu1) then var2:=var2 union {op(1,mu1)}: fi: od: if nops(var2)<>1 then F:=coeff(F,var2[1],1): fi: flo:=degree(F,P): pip:=coeff(F,P,flo): flo:=degree(pip,x): pip:=coeff(pip,x,flo): F:=F/pip: F:=numer(normal(F)): vu:=add(gu[i+1]*x^i,i=0..nops(gu)-1): halev:=add(factor(coeff(F,P,i))*P^i,i=0..degree(F,P)): vu1:=taylor(subs(P=vu,halev),x=0,nops(gu)+1): if {seq(expand(coeff(vu1,x,i)),i=1..nops(gu)-1)}<>{0} then RETURN(FAIL): fi: halev: end: Empir1:=proc(gu,x,P) local degx,degP,L,lu,che,i: for L from 1 to (nops(gu)-3)/3 do for degP from 1 to L do for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do lu:=empir(gu,degx,degP,x,P): if lu<>FAIL then che:=expand(taylor(subs(P=add(gu[i+1]*x^i,i=0..nops(gu)-1),lu),x=0,nops(gu)-1)): if {seq(coeff(che,x,i),i=0..nops(gu)-3)}<>{0} then RETURN(FAIL): else RETURN(lu): fi: fi: od: od: od: FAIL: end: with(gfun): gfun[maxdegeqn]:=14: EmpirF:=proc(gu,t,P) local eq1,eq2,vu,lu,i: if nops(gu)<20 then print(`list is too short `): RETURN(FAIL): fi: eq1:=listtoalgeq([op(1..nops(gu)-10,gu)],P(t),[ogf]): if eq1=FAIL then RETURN(FAIL): fi: eq1:=subs(P(t)=P,eq1[1]): eq2:=listtoalgeq(gu,P(t),[ogf]): if eq2=FAIL then RETURN(FAIL): fi: eq2:=subs(P(t)=P,eq2[1]): if eq1<>eq2 then RETURN(FAIL): fi: vu:=add(gu[i+1]*t^i,i=0..nops(gu)-1): lu:=taylor(subs(P=vu,eq2),t=0,nops(gu)+4): if {seq(expand(coeff(lu,t,i)),i=1..nops(gu)-2)}<>{0} then RETURN(FAIL): fi: eq2: end: Empir:=proc(gu,x,P) local i,mu: for i from 10 by 10 to 10*trunc(nops(gu)/10) do mu:=Empir1([op(1..i,gu)],x,P): if mu<>FAIL then RETURN(mu): fi: od: Empir1(gu,x,P): end: #Empirx1(gu,x,t,P,deg): inputs a list (starting at 0) of polynomials in the variable x, another variable t #tries to find an algebraic equation in t,x,P of degree deg in x #Try: #Empirx1(FirstKw({1,0,-1},0,x,20),x,t,P,3); Empirx1:=proc(gu,x,t,P,deg) local degt, degP,eq0,Eq,POL,a,b,eq,i,lu,mu: eq0:=Empir(subs(x=0,gu),t,P): if eq0=FAIL then RETURN(FAIL): else Eq[0]:=eq0: fi: eq0:=Empir(subs(x=1,gu),t,P): if eq0=FAIL then RETURN(FAIL): else Eq[1]:=eq0: fi: eq0:=Empir(subs(x=2,gu),t,P): if eq0=FAIL then RETURN(FAIL): else Eq[2]:=eq0: fi: degt:=degree(eq0,t): degP:=degree(eq0,P): for i from 3 to deg+2 do eq0:=empir(subs(x=i,gu),degt,degP,t,P): if eq0=FAIL then RETURN(FAIL): else Eq[i]:=eq0: fi: od: for a from 0 to degP do for b from 0 to degt do mu:=[seq(coeff(coeff(Eq[i],P,a),t,b),i=0..deg+2)]: lu:=GP(mu,x): if lu=FAIL then RETURN(FAIL): else POL[a,b]:=factor(lu): fi: od: od: eq:=add(add(POL[a,b]*P^a*t^b,a=0..degP),b=0..degt): eq0:=Empir(subs(x=deg+3,gu),t,P): if expand(subs(x=deg+3,eq)-eq0)<>0 then RETURN(FAIL): else eq:=add(factor(coeff(eq,P,i))*P^i,i=0..degP): RETURN(eq): fi: end: #Empirx(gu,x,t,P,MaxDeg): inputs a list (starting at 0) of polynomials in the variable x, another variable t #tries to find an algebraic equation in t,x,P of degree <=MaxDeg in x #Try: #Empirx(FirstKw({1,0,-1},0,x,20),x,t,P,3); Empirx:=proc(gu,x,t,P,MaxDeg) local deg,eq: for deg from 1 to MaxDeg do eq:=Empirx1(gu,x,t,P,deg): if eq<>FAIL then RETURN(eq): fi: od: FAIL: end: #New stuff for 2D walks #Wfe2D(S,ftxy, t,x,y,Cx,Cy,C): inputs a set of pairs of integers S, #indicating the allowed steps, # a symbol ftxy (denoting f(t,x,y), see below), variables t, and x, and y #symbols Cx,Cy, C #where Cx[i] denotes the coefficient of x^i in f(t,x,y) and #Cy[j] denotes the coefficient of y^i in f(t,x,y) and #and #C[i,j] denotes the coefficient of x^i*y^j # outputs the pair [eq,[H1,H2]], #where eq is functional equation for f(t,x,y), the generating #function of walks where the coefficient of t^i *x^j1*y^j2 is the number of ways #of walking i steps on the positive quarter-plane starting at the origin #and winding up at [j1,j2] and never going to outside the non-negative quater-plane #C[j1,j2] denotes the coefficient of x^j1*y^j2 in f(t,x,y), #H1 is the highest j1 that shows up as a C[j1,-] and H2 is the highest j2 that shows up a a C[-,j2] #Try: #Wfe2D({[1,0],[-1,0],[0,1],[0,-1]},ftxy, t,x,y,Cx,Cy,C); Wfe2D:=proc(S,ftxy,t,x,y,Cx,Cy,C) local eq,s,j1,j2,H1,H2,mu: eq:=1: H1:=-min(seq(s[1], s in S))-1: H2:=-min(seq(s[2], s in S))-1: for s in S do mu:=ftxy-add(Cx[j1]*x^j1,j1=0..-s[1]-1)-add(Cy[j2]*y^j2,j2=0..-s[2]-1)+ add(add(C[j1,j2]*x^j1*y^j2,j1=0..-s[1]-1),j2=0..-s[2]-1): eq:=eq+t*x^s[1]*y^s[2]*mu: od: [eq,[H1,H2]]: end: #NuWalks2D(Steps,t,x): Given a set of 2D steps Steps (pairs of positive or negative #integers or 0), a non-neg. integer t and a pair x returns the number of walks #(Sequences in Steps) after time t, arrive at x, and always stay in the non-negative quarter-plane #For example, try: NuWalks2D({[1,1],[-1,0],[0,-1]},12,[0,0]); NuWalks2D:=proc(Steps,t,x) local s: option remember: if t=0 then if x=[0,0] then RETURN(1): else RETURN(0): fi: fi: if x[1]<0 or x[2]<0 then RETURN(0): fi: add(NuWalks2D(Steps,t-1,[x[1]-s[1],x[2]-s[2]]), s in Steps): end: #Seq2D(Steps,N,e): the enumerating sequence for the number of walks with i steps drawn from the set of steps Steps, #starting at [0,0] and ending at the point e, and always staying in the non-negative quarter-plane #Seq2D({[1,1],[-1,-1],[1,0],[-1,0]},20,[0,0]); #Seq2D({[1,1],[-1,0],[0,-1]},21,[0,0]); Seq2D:=proc(Steps,N,e) local t: option remember: [seq(NuWalks2D(Steps,t,e),t=0..N)]: end: #OneStepFeK2D(g,F,t,x,y,ftxy,C,Cx,Cy,H,K): given a formal power series g in t, x, and y #and a functional equation F phrased in terms of C[i], i=0..H, #and a positive integer K, applies one step of the equation #g->F(g) and truncates it to K terms in t. #Try: #lu:=Wfe2D({[1,0],[-1,0],[0,1],[0,-1]},ftxy, t,x,y,Cx,Cy,C): #OneStepFeK2D(1,lu[1],t,x,y,ftxy,C,Cx,Cy,lu[2],4): OneStepFeK2D:=proc(g,F,t,x,y,ftxy,C,Cx,Cy,H,K) local gu,i,j: gu:=subs(ftxy=g,F): for i from 0 to H[1] do gu:=expand(subs(Cx[i]=coeff(g,x,i),gu)): od: for i from 0 to H[2] do gu:=expand(subs(Cy[i]=coeff(g,y,i),gu)): od: for i from 0 to H[1] do for j from 0 to H[2] do gu:=expand(subs(C[i,j]=coeff(coeff(g,x,i),y,j),gu)): od: od: add(coeff(gu,t,i)*t^i,i=0..K): end: #FirstK2D(F,t,x,y,ftxy,C,Cx,Cy,H,K): Given #a functional equation F phrased in terms of C[i,j], i=0..H[1], j=0.., H[2], and Cx[i], i=0..H[1], and Cy[i], i=0..H[2] #finds the first K terms of the formal power series in t, with polynomials in x and y, that satisfies it. Try: #lu:=Wfe2D({[1,1],[-1,0],[0,-1]},ftxy, t,x,y,Cx,Cy,C): #FirstK2D(lu[1],t,x,y,ftxy,C,Cx,Cy,lu[2],10); FirstK2D:=proc(F,t,x,y,ftxy,C,Cx,Cy,H,K) local gu,i: gu:=0: for i from 1 to K+2 do gu:=OneStepFeK2D(gu,F,t,x,y,ftxy,C,Cx,Cy,H,K): od: [seq(coeff(gu,t,i),i=0..K)]: end: #Seq2DxySlow(S,N,x,y): inputs a set of steps and finds the first N terms of the enumerating polynomials for #n-step walks in the non-negative quarter plane, according to the weight x^dest[1]*y^dest[2], where #dest is the destination. It uses the Functinal Equation, hence is slower than Seq2Dxy #Try: #Seq2DxySlow({[1,1],[-1,-1],[1,0],[-1,0]},20,x,y); #Seq2DxySlow({[1,1],[-1,0],[0,-1]},21,x,y); Seq2DxySlow:=proc(S,N,x,y) local lu,t,ftxy,C,Cx,Cy: lu:=Wfe2D(S,ftxy, t,x,y,Cx,Cy,C): FirstK2D(lu[1],t,x,y,ftxy,C,Cx,Cy,lu[2],N) : end: #Seq2Dxy(S,N,x,y): inputs a set of steps and finds the first N terms of the enumerating polynomials for #n-step walks in the non-negative quarter plane, according to the weight x^dest[1]*y^dest[2], where #dest is the destination #Try: #Seq2Dxy({[1,1],[-1,-1],[1,0],[-1,0]},20,x,y); #Seq2Dxy({[1,1],[-1,0],[0,-1]},21,x,y); Seq2Dxy:=proc(S,N,x,y) local e1,e2,n,k1,k2,i1: option remember: k1:=max(seq(S[i1][1],i1=1..nops(S))): k2:=max(seq(S[i1][2],i1=1..nops(S))): [seq(add(add(NuWalks2D(S,n,[e1,e2])*x^e1*y^e2,e1=0..k1*n),e2=0..k2*n),n=0..N)]: end: #InfoH2D(S,n,N,MaxC): #inputs a set of steps, symbols n and N (N is the shift operator in n, i.e. Na(n):=a(n+1)), and a positive #integer MaxC, outputs, if possible, the pair [Ope0,Ope1], where Ope0 is the holonomic representation #iwith ORDER+Degree<=MaxC of the sequence enumerating walks that start at [0,0], use the #step-set S, and return to [0,0], all the while staying in the first quadrant, and #Ope1 is the holonomic representation of such walks that end anywhere you want (in the first quardrant, of course). #For example, For the Kreweras walks try: #InfoH2D({[1,1],[-1,0],[0,-1]},n,N,10); #For the Gessel walks, try: #InfoH2D({[1,1],[-1,-1],[1,0],[-1,0]},n,N,10); InfoH2D:=proc(S,n,N,MaxC) local D1,K,gu,Ope0,Ope1: option remember: K:=max(seq((2+D1)*(1+MaxC-D1)+6,D1=0..MaxC)): gu:=Seq2D(S,K,[0,0]): gu:=[op(2..nops(gu),gu)]: Ope0:=Findrec(gu,n,N,MaxC): gu:=Seq2Dxy(S,K,1,1): gu:=[op(2..nops(gu),gu)]: Ope1:=Findrec(gu,n,N,MaxC): [Ope0,Ope1]: end: #InfoA2D(S,t,P,K): #inputs a set of steps, symbols t and P (t is the variable and P is the symbol for the generating function) #an integer K, outputs, if possible, the pair [POL0,POL1], where POL0 is the algebraic equation #satisfied by the generating function of sequence enumerating walks that start at [0,0], use the #step-set S, and return to [0,0], all the while staying in the first quadrant, and #POL1 is the algebaric equation of such walks that end anywhere you want (in the first quardrant, of course). #It is guessed using the first K values as "training set", If none is found it returns FAIL. #(it does not mean that there is no such equation, only that K terms do not suffice to guess it) #For example, For the Kreweras walks try: #InfoA2D({[1,1],[-1,0],[0,-1]},t,P,60); #For the Gessel walks, try: #InfoA2D({[1,1],[-1,-1],[1,0],[-1,0]},t,P,60); InfoA2D:=proc(S,t,P,K) local gu,P0,P1: gu:=Seq2D(S,K,[0,0]): P0:=Empir(gu,t,P): gu:=Seq2Dxy(S,K,1,1): P1:=Empir(gu,t,P): [P0,P1]: end: #PrintFE2D(eq,ftxy,t,x,y,C,Cx,Cy,P,Coeff):prints nicely the functional equation eq phrased in terms of ftx,t,x,C. #Try: #PrintFE2D(Wfe2D(S,ftxy, t,x,y,Cx,Cy,C),t,x,y,C,Cx,Cy,P,Coeff): PrintFE2D:=proc(eq,ftxy,t,x,y,C,Cx,Cy,P,Coeff) local H,eq1,i1,i2: H:=eq[2]: eq1:=eq[1]: eq1:=subs(ftxy=P(t,x,y),eq1): for i1 from 0 to H[1] do for i2 from 0 to H[2] do eq1:=subs(C[i1,i2]=Coeff([x,y],[i1,i2],P(t,x,y)),eq1): od: od: for i1 from 0 to H[1] do eq1:=subs(Cx[i1]=Coeff(x,i1,P(t,x,y)),eq1): od: for i2 from 0 to H[2] do eq1:=subs(Cy[i2]=Coeff(y,i2,P(t,x,y)),eq1): od: print(P(t,x,y)=eq1): end: #InfoH2Dv(S,n,MaxC): Verbose version of InfoH2D(S,n,N,MaxC) (q.v.) #For example, For the Kreweras walks try: #InfoH2Dv({[1,1],[-1,0],[0,-1]},n,10); #For the Gessel walks, try: #InfoH2Dv({[1,1],[-1,-1],[1,0],[-1,0]},n,10); InfoH2Dv:=proc(S,n,MaxC) local N,gu,f,g,ope,Ini,i,t0: t0:=time(): gu:=InfoH2D(S,n,N,MaxC): if gu[1]<>FAIL then ope:=gu[1][1]: Ini:=gu[1][2]: print(``): print(`Theorem: Let f(n) be the number of n-step walks, starting and ending at [0,0], always staying in the first quarant, i.e. the region`): print(` x>=0, y>=0, and always using one of the following steps`): print(``): print(S): print(``): print(`Then f(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0 ): print(``): print(`subject to the initial conditions`): print(``): print(seq(f(i)=Ini[i],i=1..nops(Ini))): print(``): print(`the recurrence in Maple input notation is`): print(``): lprint(add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0 ): print(``): print(`For the sake of the OEIS, here are the first 31 terms`): print(``): lprint(Seq2D(S,30,[0,0]) ): print(``): fi: if gu[2]<>FAIL then ope:=gu[2][1]: Ini:=gu[2][2]: print(`We also have the following theorem.`): print(`Theorem: Let g(n) be the number of n-step walks, starting at [0,0], always staying in the first quarant, i.e. the region`): print(` x>=0, y>=0, and ending whenever they want there, always using one of the following steps`): print(``): print(S): print(``): print(`Then g(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*g(n+i),i=0..degree(ope,N))=0 ): print(``): print(`subject to the initial conditions`): print(``): print(seq(g(i)=Ini[i],i=1..nops(Ini))): print(``): print(`the recurrence in Maple input notation is`): print(``): lprint(add(coeff(ope,N,i)*g(n+i),i=0..degree(ope,N))=0 ): print(``): print(`For the sake of the OEIS, here are the first 31 terms`): print(``): lprint(Seq2Dxy(S,30,1,1) ): print(``): fi: print(`This took`, time()-t0, `seconds. `): end: #Tikva(S,MaxC): inputs a set of 2-dimensional steps, S, a psoitive integer MaxC, and a positive integer L0 #finds out whether the generating function of the first-quadrant walks f_S(x,y) is such that #f_S(x,0) and f_S(0,y) are holonomic in x,y. It returns #Try: #Tikva({[1,0],[0,1],[-1,-1]},10); Tikva:=proc(S,MaxC) local D1, K, n,N,gu,ope,lu,s: option remember: if min(seq(op(s),s in S))<-1 then print(`Not yet implented`): RETURN(FAIL): fi: K:=max(seq((2+D1)*(1+MaxC-D1)+6,D1=0..MaxC)): gu:=Seq2Dxy(S,K,11/3,0): ope:=Findrec(gu,n,N,MaxC): lu:=[]: if ope<>FAIL then lu:=[op(lu), [degree(numer(ope[1]),n),degree(ope[1],N)]]: else lu:=[op(lu),FAIL]: fi: gu:=Seq2Dxy(S,K,0,14/5): ope:=Findrec(gu,n,N,MaxC): if ope<>FAIL then lu:=[op(lu), [degree(numer(ope[1]),n),degree(ope[1],N)]]: else lu:=[op(lu),FAIL]: fi: lu: end: #InfoH2Dv1(S,n,MaxC,n0): Verbose version of InfoH2D(S,n,N,MaxC) (q.v.) #but also makes it a numbered Theorem #For example, For the Kreweras walks try: #InfoH2Dv1({[1,1],[-1,0],[0,-1]},n,10,n0); #For the Gessel walks, try: #InfoH2Dv1({[1,1],[-1,-1],[1,0],[-1,0]},n,10,5); InfoH2Dv1:=proc(S,n,MaxC,n0) local N,gu,f,g,ope,Ini,i: gu:=InfoH2D(S,n,N,MaxC): if gu[1]<>FAIL then ope:=gu[1][1]: Ini:=gu[1][2]: print(``): print(`Theorem `, n0, `:Let f(n) be the number of n-step walks, starting and ending at [0,0], always staying in the first quarant, i.e. the region`): print(` x>=0, y>=0, and always using one of the following steps`): print(``): print(S): print(``): print(`Then f(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0 ): print(``): print(`subject to the initial conditions`): print(``): print(seq(f(i)=Ini[i],i=1..nops(Ini))): print(``): print(`the recurrence in Maple input notation is`): print(``): lprint(add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0 ): print(``): print(`For the sake of the OEIS, here are the first 31 terms`): print(``): lprint(Seq2D(S,30,[0,0]) ): print(``): fi: if gu[2]<>FAIL then ope:=gu[2][1]: Ini:=gu[2][2]: print(`We also have the following theorem.`): print(`Theorem `, cat(n0,A), `: Let g(n) be the number of n-step walks, starting at [0,0], always staying in the first quarant, i.e. the region`): print(` x>=0, y>=0, and ending whenever they want there, and always using one of the following steps`): print(``): print(S): print(``): print(`Then g(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*g(n+i),i=0..degree(ope,N))=0 ): print(``): print(`subject to the initial conditions`): print(``): print(seq(g(i)=Ini[i],i=1..nops(Ini))): print(``): print(`the recurrence in Maple input notation is`): print(``): lprint(add(coeff(ope,N,i)*g(n+i),i=0..degree(ope,N))=0 ): print(``): print(`For the sake of the OEIS, here are the first 31 terms`): print(``): lprint(Seq2Dxy(S,30,1,1) ): print(``): fi: end: #Sefer(k,MaxC): inputs a positive integer k between 3 and 8 and outputs all the theorems #about walks with k small steps, where the steps are subsets of #{[1,0],[-1,0],[0,1],[0,-1],[1,1],[1,-1],[-1,1],[-1,-1]}. Try #Sefer(3,10); Sefer:=proc(k,MaxC) local G,lu,S,t0,c1,c2,efes,n0,vu, N, n,gu,tov1,tov2, akshan,vu1: t0:=time(): efes:=0: c1:=0: c2:=0: n0:=0: tov1:=0: tov2:=0: akshan:={}: G:={[1,0],[-1,0],[0,1],[0,-1],[1,1],[1,-1],[-1,1],[-1,-1]}: print(``): print(`Systematic Enumeration of Two-Dimensional Lattice Walks Confined to the First Quadrant with`, k, ` Short Steps`): print(`With Holonomic Complexity Not Exceeding`, MaxC): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`In this article we will attempt to enumerate, in the context of the Holonomic Ansatz,`): print(`lattice walks in the Two Dimensional Non-negative Square Lattice N^2 that start at the origin, [0,0],`): print(`(in other words an infinite chessboard, starting at the bottom-left corner)`): print(`where N={0,1,2, ...} is the set of non-negative integers. `): print(`We will try to do it for all possible subsets with`, k, ` elements of the 2^8=256 set of king-steps {[+-1,+-1], and [1,0],[0,1],[0,1],[1,0]},`): print(`and try and find holonomic representations, i.e., linear recurrence equations with polynomial coefficients whose complexity is <=`, MaxC, `,`): print(`in other words, the order of the recurrence plus the highest degree of the coefficients of the recurrence, in the discrete variable, n, is <=`, MaxC): print(`We will try to enumerate both those walks that start and return to the origin, and the total number of walks,`): print(` (still starting at the origin) regardless of final destination.`): print(`We will discard obviously trivial cases, where the enumerating sequence is eventually a constant.`): print(`This work reproduces, from the perspective of Doron Zeilberger's semi-rigorous and non-rigorous philosophy,`): print(` interesting, "rigorous" investigations by smart humans like Marni Mishna, Mireille Bousquet-Melou, Manuel Kauers, Christoph Koutschan`): print(`my beloved master Dr. Z. , and other people. `): lu:=choose(G,k): print(``): print(`-------------------------------------------------------------------------`): print(``): efes:=0: c1:=0: c2:=0: n0:=0: for S in lu do print(`Invesigating the set of steps`, S): vu1:= Seq2Dxy(S,15,1,1): vu:= Seq2D(S,15,[0,0]): if nops({op(6..nops(vu1),vu1)})=1 or nops({op(6..nops(vu),vu)})=1 then efes:=efes+1: print(`This case is trivial`): else gu:=InfoH2D(S,n,N,MaxC): if gu[1]=FAIL then c1:=c1+1: else tov1:=tov1+1: fi: if gu[2]=FAIL then akshan:=akshan union {S}: c2:=c2+1: else tov2:=tov2+1: fi: if gu[1]<>FAIL or gu[2]<>FAIL then n0:=n0+1: print(`We have `): InfoH2Dv1(S,n,MaxC,n0): fi: fi: od: print(``): print(`----------------------------------------------`): print(``): print(`We have discovered`, tov1+tov2, `exciting theorems in enumerative combinatorics. We didn't bother with proofs, since we know how it could be done`): print(` so why bother?`): print(``): print(`We note that`, efes , `cases were trivial`, c1, `cases do not have a holonomic representation of complexity <=`, MaxC, `for the enuemrating sequence`): print(`of walks that start and end at the origin, and`, c2, `cases do not have such representations for enumerating walks that start at the origin`): print(`and end anywhere in the first quardrat (of course staying all the time in that quardrant.)`): print(``): print(`We also note that the enumerating sequence for walks that end anywhere with the following set of set of steps `): lprint(akshan): print(`do not have a holonomic representation of complexity <=`, MaxC, ` and are possibly non-holonomic. `): print(``): print(`----------------------------------------------`): print(`The whole thing took`, time()-t0, `seconds. `): end: #PrintFE2Ds(eq,ftxy,t,x,y,C,Cx,Cy,P,Coeff):prints nicely the functional equation eq phrased in terms of ftx,t,x,C. #Try: #PrintFE2Ds(Wfe2D(S,ftxy, t,x,y,Cx,Cy,C),t,x,y,C,Cx,Cy,P,Coeff): PrintFE2Ds:=proc(eq,ftxy,t,x,y,C,Cx,Cy,P,Coeff) local H,eq1,i1,i2: H:=eq[2]: eq1:=eq[1]: eq1:=subs(ftxy=P(t,x,y),eq1): for i1 from 0 to H[1] do for i2 from 0 to H[2] do if i1=0 and i2=0 then eq1:=subs(C[i1,i2]=P(t,0,0),eq1): else eq1:=subs(C[i1,i2]=Coeff([x,y],[i1,i2],P(t,x,y)),eq1): fi: od: od: for i1 from 0 to H[1] do if i1=0 then eq1:=subs(Cx[i1]=P(t,0,y),eq1): else eq1:=subs(Cx[i1]=Coeff(x,i1,P(t,x,y)),eq1): fi: od: for i2 from 0 to H[2] do if i2=0 then eq1:=subs(Cy[i2]=P(t,x,0),eq1): else eq1:=subs(Cy[i2]=Coeff(y,i2,P(t,x,y)),eq1): fi: od: print(P(t,x,y)=eq1): end: #InfoH2DvWP(S,n,MaxC): Verbose version of InfoH2D(S,n,N,MaxC) (q.v.) with a sketch of proof #For example, For the Kreweras walks try: #InfoH2DvWP({[1,1],[-1,0],[0,-1]},n,10); InfoH2DvWP:=proc(S,n,MaxC) local N,gu,f,g,ope,Ini,i,t0,lu,w,a,b,F,ka,ftxy, t,x,y,Cx,Cy,C,eq,ka1: t0:=time(): gu:=InfoH2D(S,n,N,MaxC): if gu[1]=FAIL or gu[2]=FAIL then print(`Could not find recurrences of complexity <=`, MaxC, `for BOTH enumerating sequences for walks that start and end at the origin`): print(`and walks that start at the origin and end everywhere with set of steps`, S): RETURN(FAIL): fi: lu:=Tikva(S,MaxC): if lu[1]=FAIL or lu[2]=FAIL then print(`Could find theorems but could not find proofs, but I am sure that there is one, just make MaxC bigger than`, MaxC, `and try again. `): RETURN(FAIL): fi: ope:=gu[1][1]: Ini:=gu[1][2]: print(``): print(`Theorem: Let f(n) be the number of n-step walks, starting and ending at [0,0], always staying in the first quarant, i.e. the region`): print(` x>=0, y>=0, and always using one of the following steps`): print(``): print(S): print(``): print(`Then f(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0 ): print(``): print(`subject to the initial conditions`): print(``): print(seq(f(i)=Ini[i],i=1..nops(Ini))): print(``): print(`the recurrence in Maple input notation is`): print(``): lprint(add(coeff(ope,N,i)*f(n+i),i=0..degree(ope,N))=0 ): print(``): print(`For the sake of the OEIS, here are the first 31 terms`): print(``): lprint(Seq2D(S,30,[0,0]) ): print(``): ope:=gu[2][1]: Ini:=gu[2][2]: print(`We also have the following theorem.`): print(`Theorem: Let g(n) be the number of n-step walks, starting at [0,0], always staying in the first quarant, i.e. the region`): print(` x>=0, y>=0, and ending whenever they want there, always using one of the following steps`): print(``): print(S): print(``): print(`Then g(n) satisfies the following linear recurrence equation with polynomial coefficients`): print(``): print(add(coeff(ope,N,i)*g(n+i),i=0..degree(ope,N))=0 ): print(``): print(`subject to the initial conditions`): print(``): print(seq(g(i)=Ini[i],i=1..nops(Ini))): print(``): print(`the recurrence in Maple input notation is`): print(``): lprint(add(coeff(ope,N,i)*g(n+i),i=0..degree(ope,N))=0 ): print(``): print(`For the sake of the OEIS, here are the first 31 terms`): print(``): lprint(Seq2Dxy(S,30,1,1) ): print(``): print(`----------------------------------`): print(``): print(`Sketch of proof`): print(``): print(`Let w(n,a,b) be the number of walks of length n with step-set`, S, `that start at [0,0] and end at [a,b] `): print(`and always staying in the first-quadrant. Let F(t,x,y) be the three-variable generating function, i.e.:`): print(``): print(F(t,x,y)=Sum(Sum(Sum(w(n,a,b)*t^n*x^a*y^b,b=0..infinity),a=0..infinity),n=0..infinity)): print(``): print(`It is readily seen that F(t,x,y) satisfies the following Functional Equation that uniquely determines it: `): eq:=Wfe2D(S,ftxy, t,x,y,Cx,Cy,C): PrintFE2Ds(eq,ftxy,t,x,y,C,Cx,Cy,F,Coeff): eq:=ftxy-eq[1]: ka:=solve(eq,ftxy): ka1:=subs({Cx[0]=F(t,x,0),Cy[0]=F(t,0,y)},ka): print(`Solving for F(t,x,y), this is equivalent to`): print(F(t,x,y)=ka1): print(`It is possible to guess, from the initial values, a linear recurrence equation with polynomial coefficients (in t and x) `): print(`for the sequence of polynomials in x given by the coefficients, in t, of F(t,x,0)`): print(`of order`, lu[1][1], `and degree of coefficients`, lu[1][2] ): print(`This immediately implies that F(t,x,0) is most probably D-finite and it is easy to find the linear differential equation in t, with polynomial coefficients`): print(`in t and x (so far conjecturally) satisfied by it. `): print(`It is also possible to guess, from the initial values, a linear recurrence equation with polynomial coefficients (in t and y) `): print(`for the sequence of polynomials in y given by the coefficients, in t, of F(t,0,y)`): print(`of order`, lu[2][1], `and degree of coefficients`, lu[2][2] ): print(`This immediately implies that F(t,0,y) is most probably D-finite and it is easy to find the linear differential equation in t, with polynomial coefficients`): print(`in t and y (so far conjecturally) satisfied by it. `): print(`Now let's go backwards! Let f1(t,x) and f2(t,y) be the unique formal power series in t that satisfy the above-mentioned guessed differential`): print(`equation and define `): print(``): print(G(t,x,y)=subs({Cx[0]=f1(t,x),Cy[0]=f2(t,y)},ka)): print(`it is purely routine to find a differential equation with polynomial coefficients satisfied by G(t,x,y).`): print(`It is also purely routine to check (still in the holonomic ansatz in the single-variable t,`): print(`that the defining functional equation for F(t,x,y) is satisfied by the newly arrived G(t,x,y).`): print(`Hence by uniqueness, F(t,x,y)=G(t,x,y), and the above-mentioned differential equation in t is satisfied by F(t,x,y).`): print(`Now plugging-in x=0, y=0, gives us a differential equation satisfied by F(t,0,0) and plugging-in x=1, y=1 gives us a differential equation`): print(`satisfied by F(t,1,1) that translate to the recurrenes claimed in the above two theorems. QED! (modulo purely routine calculations)`): print(`This took`, time()-t0, `seconds. `): end: #TikvaLight(S,MaxC): inputs a set of 2-dimensional steps, S, a psoitive integer MaxC, and a positive integer L0 #finds out whether the generating function of the first-quadrant walks f_S(x,y) is such that #f_S(x,0) and f_S(0,y) are holonomic in x,y, but by plugging in x=2 It returns #Try: #TikvaLight({[1,0],[0,1],[-1,-1]},10); TikvaLight:=proc(S,MaxC) local D1, K, n,N,gu,ope,s: option remember: if min(seq(op(s),s in S))<-1 then print(`Not yet implented`): RETURN(FAIL): fi: K:=max(seq((2+D1)*(1+MaxC-D1)+6,D1=0..MaxC)): gu:=Seq2Dxy(S,K,2,0): ope:=Findrec(gu,n,N,MaxC): if ope<>FAIL then [degree(numer(ope[1]),n),degree(ope[1],N)]: else FAIL: fi: end: