print(`This is the package POLYMER`): print(`written by Arvind Ayyer and Doron Zeilberger`): print(`It is written for the analysis of directed lattice walks`): print(`in two dimensions with arbitrary set of steps`): print(`It uses some functions from GuessHolo2`): print(`written by D. Zeilberger`): print(`To see the list of procedures, type Help();`): print(`For help on a specific procedure, type Help(procedure_name);`): Help := proc() if args=NULL then print(`Contains the following procedures`): print(`Number of walks: polymer, WEpolymer, polymerBE, WEpolymerBE`): print(`Recurrences: stepsrec, WEstepsrec,`): print(`Generating functions: rigorgf, rigorgfWE, RGF2D, RGF2DWE`): print(`Plotting commands: plotFE, plotFEWE, ForceWE, plotinfFEWE`): print(`NOTE: The plots are best viewed in the Graphical Version`): fi: if nops([args])=1 and op(1,[args])=`polymer` then print(`polymer(Steps,a,w) uses steps given in Steps `): print(`to calculate the number of ways of walking upto a from the origin`): print(`where the walk is confined between the lines x-y=0 and x-y=w.`): print(`For example, try polymer({[0,1],[1,0]},[2,2],3);`): fi: if nops([args])=1 and op(1,[args])=`polymerBE` then print(`polymer(Steps,a,b,w) uses steps given in Steps `): print(`to calculate the number of ways of walking from a to b`): print(`where the walk is confined between the lines x-y=0 and x-y=w.`): print(`For example, try polymerBE({[0,1],[1,0]},[1,1],[2,2],3);`): fi: if nops([args])=1 and op(1,[args])=`WEpolymer` then print(`WEpolymer(Steps,a,w,t,s) uses steps given in Steps `): print(`to calculate the enumerated walk from 0 to a`): print(`where the walk is confined between the lines x-y=0 and x-y=w`): print(`and t,s mark the number of times they are hit respectively.`): print(`For example, try WEpolymer({[0,1],[1,0]},[2,2],2,t,s);`): fi: if nops([args])=1 and op(1,[args])=`WEpolymerBE` then print(`WEpolymer(Steps,a,b,w,t,s) uses steps given in Steps `): print(`to calculate the enumerated walk from a to b`): print(`where the walk is confined between the lines x-y=0 and x-y=w`): print(`and t,s mark the number of times they are hit respectively.`): print(`For example, try WEpolymerBE({[0,1],[1,0]},[1,1],[2,2],2,t,s);`): fi: if nops([args])=1 and op(1,[args])=`stepsrec` then print(`stepsrec(Steps,METADEGREE,METAORDER,N,w,W) finds an operator, if it exists,`): print(`of order METAORDER and degree METADEGREE which annihilates denominators`): print(`of subsequent generating functions as a function of N of sufficiently many widths.`): print(`For example, try stepsrec({[1,0],[0,1]},0,2,N,w,W);`): fi: if nops([args])=1 and op(1,[args])=`WEstepsrec` then print(`WEstepsrec(Steps,METADEGREE,METAORDER,N,w,W) finds an operator, if it exists,`): print(`of order METAORDER and degree METADEGREE which annihilates denominators`): print(`of subsequent weight enumerating generating functions as a function of N,t,s of sufficiently many widths.`): print(`For example, try WEstepsrec({[1,0],[0,1]},0,2,N,w,W,t,s);`): fi: if nops([args])=1 and op(1,[args])=`rigorgf` then print(`rigorgf(Steps,w,z) calculates the generating function exactly`): print(`of walks with Steps and width w in the variable z.`): print(`Try, for example, rigorgf({[0,1],[1,0]},2,z);`): fi: if nops([args])=1 and op(1,[args])=`rigorgfWE` then print(`rigorgfWE(Steps,w,z,t,s) calculates the generating function exactly`): print(`of walks with Steps and width w in the variable z`): print(`where t,s mark the number of times the walk hits`): print(`the lines x-y=0 and x-y=w respectively.`): print(`Try, for example, rigorgfWE({[0,1],[1,0]},2,z,t,s);`): fi: if nops([args])=1 and op(1,[args])=`RGF2D` then print(`RGF2D(Steps,z,F) rigorously finds the polynomial equation`): print(`satisfied by the generating function F(z)`): print(`for a walk with arbitrary sequence of steps with infinite width.`): print(`Try, for example, RGF2D({[0,1],[1,0]},z,F);`): fi: if nops([args])=1 and op(1,[args])=`RGF2DWE` then print(`RGF2DWE(Steps,z,F,t) rigorously finds the polynomial equation`): print(`satisfied by the generating function F(z;t)`): print(`for a walk with arbitrary sequence of steps with infinite width`): print(`where the parameter t marks the number of times the walk`): print(`hits the line x-y=0.`): print(`Try, for example, RGF2DWE({[0,1],[1,0]},z,F,t);`): fi: if nops([args])=1 and op(1,[args])=`rgfprove` then print(`rgfprove(Steps,OPE,W,N) when given OPE(W,N), the ope satisfied by the denominators`): print(`of the generating functions as a function of N and the shift operator W`): print(`verifies whether the implied relation for the GF is true.`): print(`Try, for example, rgfprove({[0,1],[1,0]},N-W+W^2,W,N)`): fi: if nops([args])=1 and op(1,[args])=`plotFE` then print(`plotFE(Steps,winit,wfin) plots the numerical value of the free energy`): print(`when given the set of steps and the initial and final integer widths.`): print(`Try, for example, plotFE({[0,1],[1,0]},1,10); `): fi: if nops([args])=1 and op(1,[args])=`plotFEWE` then print(`plotFEWE(Steps,w,tmin,tmax,smin,smax) plots the free energy for a`): print(`given set of steps and width. The ranges of the weight enumerating`): print(`parameters t,s are to be specified.`): print(`Try, for example, plotFEWE({[0,1],[1,0]},3,1,4,1,4); `): fi: if nops([args])=1 and op(1,[args])=`ForceWE` then print(`ForceWE(Steps,w,tmin,tmax,smin,smax) plot the force for a`): print(`given set of steps and width. The ranges of the weight enumerating`): print(`parameters t,s are to be specified.`): print(`Try, for example, ForceWE({[0,1],[1,0]},3,1,4,1,4); `): fi: if nops([args])=1 and op(1,[args])=`plotinfFEWE` then print(`plotinfFEWE(Steps,tmin,tmax) finds the Free Energy for the`): print(`set Steps with the weight enumerating parameter varying from `): print(`integers tmin to tmax.`): print(`Try, for example, plotinfFEWE({[0,1],[1,0]},1,5);`): fi: end: with(combinat): with(gfun): with(SolveTools): with(plots): #polymer(Steps,a,w) uses steps given in Steps to calculate the number of ways of walking upto a where the walk is confined between the lines x-y=0 and x-y=w polymer := proc(Steps, a, w) local i,j,k,Prev: option remember: k := nops(a): if k <> nops(Steps[1]) then return FAIL: fi: if a=[0$k] then return 1: fi: #HAVE TO REPLACE THIS WITH THE PROPER CONDITION #if min(op(a))<0 then # return 0: #fi: #THIS SEEMS RIGHT if sum(a[i],i=1..nops(a)) < 0 then return 0: fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then return 0: fi: if max(seq(a[i]-a[i+1],i=1..nops(a)-1))>w then return 0: fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return add(polymer(Steps,Prev[i],w),i=1..nops(Prev)): end: #polymerBE(Steps,a,b,w) uses steps given in Steps to calculate the number of ways of walking from a to b where the walk is confined between the lines x-y=0 and x-y=w polymerBE := proc(Steps, a, b, w) local i,j,k,ba,Prev: option remember: k := nops(a): ba := b-a: if k <> nops(Steps[1]) then return FAIL: fi: #HAVE TO REPLACE THIS WITH THE PROPER CONDITION #if min(op(ba))<0 then # return 0: #fi: #THIS SEEMS RIGHT if sum(ba[i],i=1..nops(ba)) < 0 then return 0: fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then return 0: fi: if max(seq(a[i]-a[i+1],i=1..nops(a)-1))>w then return 0: fi: if min(seq(b[i]-b[i+1],i=1..nops(b)-1))<0 then return 0: fi: if max(seq(b[i]-b[i+1],i=1..nops(b)-1))>w then return 0: fi: if ba=[0$k] then return 1: fi: Prev:={seq([seq(b[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return add(polymerBE(Steps,a,Prev[i],w),i=1..nops(Prev)): end: #WEpolymerBE(Steps,a,b,w,t,s) uses steps given in Steps to calculate the number of ways of walking from a to b where the walk is confined between the lines x-y=0 and x-y=w WEpolymerBE := proc(Steps, a, b, w,t,s) local i,j,k,ba,Prev: option remember: k := nops(a): ba := b-a: if k <> nops(Steps[1]) then return FAIL: fi: if min(op(ba))<0 then return 0: fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then return 0: fi: if max(seq(a[i]-a[i+1],i=1..nops(a)-1))>w then return 0: fi: if min(seq(b[i]-b[i+1],i=1..nops(b)-1))<0 then return 0: fi: if max(seq(b[i]-b[i+1],i=1..nops(b)-1))>w then return 0: fi: if ba=[0$k] then return 1: fi: if min(seq(b[i]-b[i+1],i=1..nops(b)-1))=0 then Prev:={seq([seq(b[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return simplify(t*add(WEpolymerBE(Steps,a,Prev[i],w,t,s),i=1..nops(Prev))): fi: if max(seq(b[i]-b[i+1],i=1..nops(b)-1))=w then Prev:={seq([seq(b[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return simplify(s*add(WEpolymerBE(Steps,a,Prev[i],w,t,s),i=1..nops(Prev))): fi: Prev:={seq([seq(b[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return add(WEpolymerBE(Steps,a,Prev[i],w,t,s),i=1..nops(Prev)): end: #WEpolymer(Steps,a,w,t,s) uses steps given in Steps to calculate the weight enumerator of the number of ways of walking upto a where the walk is confined between the lines x-y=0 and x-y=w (t marks the 0-hyperplane and s marks the w-hyperplane) WEpolymer := proc(Steps, a, w,t,s) local i,j,k,Prev: option remember: k := nops(a): if k <> nops(Steps[1]) then return FAIL: fi: 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: if max(seq(a[i]-a[i+1],i=1..nops(a)-1))>w then return 0: fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))=0 then Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return simplify(t*add(WEpolymer(Steps,Prev[i],w,t,s),i=1..nops(Prev))): fi: if max(seq(a[i]-a[i+1],i=1..nops(a)-1))=w then Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return simplify(s*add(WEpolymer(Steps,Prev[i],w,t,s),i=1..nops(Prev))): fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: return add(WEpolymer(Steps,Prev[i],w,t,s),i=1..nops(Prev)): end: #stepsrec(Steps,METADEGREE,METAORDER,N,w,W) finds an operator, if it exists, of order METAORDER and degree METADEGREE which annihilates denominators of subsequent generating functions (functions of N) of sufficiently many widths stepsrec := proc(Steps,METADEGREE,METAORDER,N,w,W) local i,ans,sol,width,make_sol: width := (1+METADEGREE)*(1+METAORDER)+5+METAORDER: make_sol := proc(i) local res; res := denom(rigorgf(Steps,i,N)); `if`(coeff(res, N, 0)<0, -1, 1)*res; end; sol := Vector(width, make_sol); ans := findrec(sol,METADEGREE,METAORDER,w,W): #print(sol): RETURN(`Meta-recursion yields`,ans): end: #WEstepsrec(Steps,METADEGREE,METAORDER,N,w,W,t,s) finds an operator, if it exists, of order METAORDER and degree METADEGREE which annihilates denominators of subsequent weight enumerators (functions of N,t,s) of sufficiently many widths WEstepsrec := proc(Steps,METADEGREE,METAORDER,N,w,W,t,s) local i,ans,sol,width: width := (1+METADEGREE)*(1+METAORDER)+5+METAORDER: sol := [seq(denom(rigorgfWE(Steps,i,N,t,s)),i=1..width)]: for i from 1 to nops(sol) do if coeff(sol[i],N,0) < 0 then sol[i] := -sol[i]: fi: od: sol := Vector(sol): ans := findrec(sol,METADEGREE,METAORDER,w,W): #print(sol): RETURN(`Meta-recursion yields`,ans): end: #algeqtorec(eqn,F,z,N,n) returns the recurrence operator as a function of N and n given the algebraic equation eqn(z,F(z))=0 algeqtorec := proc(eqn,F,z,N,n) local i,j,deqns,deqn,reqns,reqn,ord,ope,P: deqns := algeqtodiffeq(eqn,F(z)): if type(deqns,set) then for i from 1 to nops(deqns) do if type(deqns[i],`+`) then deqn := deqns[i]: fi: od: else deqn := deqns: fi: reqns := diffeqtorec(deqn,F(z),P(n)): if type(reqns,set) then for i from 1 to nops(reqns) do if type(reqns[i],`+`) then reqn := reqns[i]: fi: od: else reqn := reqns: fi: #ord := nops(reqn)-1: ord := max(seq(op(1,op(2,op(i,reqn))),i=1..nops(reqn)))-n: ope := subs({seq(P(n+i)=N^i,i=0..ord)},reqn): return ope: end: #Asy(eq,F,z) given the algebraic equation satisfied by the generating function, returns the free energy Asy := proc(eqn,F,z) local i,deqn,deq,deq2,poly,sols: deqn := algeqtodiffeq(eqn,F(z)): if type(deqn,set) then for i from 1 to nops(deqn) do if type(deqn[i],`+`) then deq := deqn[i]: fi: od: else deq := deqn: fi: deq := subs({seq(diff(F(z),z$i)=D^i,i=1..nops(deq))},deq): poly := coeff(deq,D,degree(deq,D)): sols := {solve(poly,z)} minus {0}: return 1/min(seq(abs(sols[i]),i=1..nops(sols))): end: #FE(ope) given the denominator of the rational generating function calculates the free energy FE := proc(ope) local i,approx,z,soln: approx := solve(ope): if nops({approx})>1 then soln := seq(abs(Re(evalf(approx[i]))),i=1..nops({approx})): else soln := approx: fi: return log(1/min(soln)): end: #plotFE(Steps,winit,wfin) plots the free energies from w=winit to w=wfin using listplot plotFE := proc(Steps,winit,wfin) local i,llist,w,ope,asymp,F,z: llist := [seq([w,FE(denom(rigorgf(Steps,w,z)))],w=winit..wfin)]: print(llist): ope := RGF2D(Steps,z,F): asymp := Asy(ope,F,z): print(`The asymptotic value is`,evalf(log(asymp))): listplot(llist,labels=["Width","Free Energy"]): end: #plotFEWE(Steps,w,tmin,tmax,smin,smax) given the steps Steps, width w and weight enumerators t,s calculates the free energy using Asy in GuessHolo2 plotFEWE := proc(Steps,w,tmin,tmax,smin,smax) local i,j,walks,ope,approx,z,Z,t,s: approx := [seq([seq(FE(denom(rigorgfWE(Steps,w,z,i,j))),j=smin..smax)],i=tmin..tmax)]: #return approx: listplot3d(approx,labels=["t","s","Free Energy"],axes=boxed,orientation=[-74,73]): end: #ForceWE(Steps,w,tmin,tmax,smin,smax) calculates the numerical value of the force for width w and values of t and s ForceWE := proc(Steps,w,tmin,tmax,smin,smax) local t,s,ope1,ope2,z,Z,approx1,approx2: ope1 := rigorgfWE(Steps,w,z,j,i): ope2 := rigorgfWE(Steps,w+1,z,j,i): approx1 := [seq([seq(FE(subs({j=t,i=s},denom(ope1))),s=smin..smax)],t=tmin..tmax)]: approx2 := [seq([seq(FE(subs({j=t,i=s},denom(ope2))),s=smin..smax)],t=tmin..tmax)]: listplot3d(approx2-approx1,labels=["t","s","Force"],axes=boxed,orientation=[-74,73]): end: #plotinfFEWE(Steps,tmin,tmax) plots the free energy in the infinite width case in the given range of t. plotinfFEWE := proc(Steps,tmin,tmax) local llist,i,j,t,s,n,w: #Far off point from the origin n := 100: #High enough width to get reasonable answers w := 1000: llist := [seq([t,evalf(log(WEpolymer(Steps,[n+1,n+1],w,t,s)/WEpolymer(Steps,[n,n],w,t,s)))],t=tmin..tmax)]: #eqn := RGF2DWE(Steps,z,H,t): #llist := [seq([t,evalf(log(Asy(RGF2DWE(Steps,z,H,t),H,z)))],t=tmin..tmax)]: listplot(llist,labels=["t","Free Energy"]): end: #anmat(M,DEGREE,ORDER,n,N) given a list of lists M tries to find a holonomic ansatz for the non-zero part of the columns, ie considering the coefficient of N^0,N^1,etc. anmat := proc(M,DEGREE,ORDER,n,N) local i,j,k,l,x,ans,sol,rl,cl,marker: rl := nops(M): cl := nops(M[1]): sol := []: #marker[i] is the location in the i'th column of the location where the first non-zero term occurs marker := []: for i from 1 to cl do for j from 1 to rl-1 do if M[1][i] <> 0 then marker := [op(marker),1]: break: fi: if M[rl-1][i] = 0 and M[rl][i] = 0 then marker := [op(marker),rl]: break: fi: if M[j][i] = 0 and M[j+1][i] <> 0 then marker := [op(marker),j+1]: break: fi: od: od: #return marker: for k from 1 to cl do x := 0: for i from 1 to ORDER do for j from 0 to DEGREE do if x = 0 and (1+DEGREE)*(1+ORDER)+5+ORDER <= rl-marker[k] then ans := findrec( [seq(M[marker[k]+l][k],l=0..rl-marker[k])], j,i,n,N): if ans <> FAIL then sol := [op(sol),ans]: x := 1: fi: if j = DEGREE and i = ORDER and ans = FAIL then print(`Need a higher DEGREE/ORDER for a solution for N^`,k-1): x := 1: fi: fi: od: od: od: return sol: end: #rigorgf(Steps,w,z) calculates the generating function of walks with Steps and width w rigorgf := proc(Steps,w,z) local i,j,temp,scri,eqs,var,sol: eqs := {}: for i from 1 to w+1 do temp := scri[i]: for j from 1 to nops(Steps) do if i+Steps[j][2]-Steps[j][1] <= w+1 and i+Steps[j][2]-Steps[j][1] >=1 then temp := temp- z^Steps[j][2]*scri[i+Steps[j][2]-Steps[j][1]]: fi: od: if i = 1 then temp := temp -1: fi: eqs := eqs union {temp}: od: var := {}: for i from 1 to w+1 do var := var union {scri[i]}: od: sol := subs(solve(eqs,var),scri[1]): #print(`The first few terms of the Taylor series of the solution`): #print(`are`, seq(coeftayl(sol,z=0,n),n=0..10)): #print(`And the first few terms are precisely`): #print(seq(polymer(Steps,[n,n],w),n=0..10)): return sol: end: #rigorgfWE(Steps,w,z,t,s) calculates the generating function of walks with Steps and width w rigorgfWE := proc(Steps,w,z,t,s) local i,j,temp,scri,eqs,var: eqs := {}: for i from 1 to w+1 do temp := scri[i]: for j from 1 to nops(Steps) do if i+Steps[j][2]-Steps[j][1] <= w+1 and i+Steps[j][2]-Steps[j][1] >=1 then if i = 1 then temp := temp- t*z^Steps[j][2]*scri[i+Steps[j][2]-Steps[j][1]]: elif i=w+1 then temp := temp- s*z^Steps[j][2]*scri[i+Steps[j][2]-Steps[j][1]]: else temp := temp- z^Steps[j][2]*scri[i+Steps[j][2]-Steps[j][1]]: fi: fi: od: if i = 1 then temp := temp -1: fi: eqs := eqs union {temp}: od: #return eqs: var := {}: for i from 1 to w+1 do var := var union {scri[i]}: od: #return solve(eqs,var): return subs(solve(eqs,var),scri[1]): end: #RGF2D(Steps,z,F) finds the rigorous generating function for an arbitrary sequence of steps with infinite width RGF2D:= proc(Steps, z, F) local i,j,k,G,vars,eqs,maxL,temp,num,sol,diffop,recop: num := nops(Steps): maxL := 0: for i from 1 to num do if abs(Steps[i][1]-Steps[i][2]) > maxL+1 then maxL := abs(Steps[i][1]-Steps[i][2])-1: fi: od: eqs := F[1][1]-1-G[1][1]*F[1][1]: for i from 1 to nops(Steps) do if Steps[i][1] = Steps[i][2] then eqs := eqs-z^Steps[i][1]*F[1][1]: fi: od: eqs := {eqs}: temp := G[1][1]: for i from 1 to num do for j from 1 to num do if Steps[i][1] > Steps[i][2] and Steps[j][1] < Steps[j][2] then temp := temp - z^(Steps[j][2]+Steps[i][2])*F[abs(Steps[i][1]-Steps[i][2])][abs(Steps[j][1]-Steps[j][2])]: fi: od: od: eqs := eqs union {temp}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do if i <> 1 and j <> 1 then if i >= j then eqs := eqs union {G[i][j]-G[i-j+1][1]}: else eqs := eqs union {G[i][j]-G[1][j-i+1]}: fi: fi: od: od: for i from 1 to maxL+1 do for j from 1 to maxL+1 do if i <> j then if i <> 1 or j <> 1 then temp := F[i][j]: for k from 1 to min(i,j) do if i < j then temp := temp - F[i-k+1][1]*G[k][j]: else temp := temp - G[i][k]*F[1][j-k+1]: fi: od: eqs := eqs union {temp}: fi: fi: od: od: for k from 2 to maxL+1 do temp := F[k][k]-F[1][1]: for i from 2 to k do temp := temp - G[k-i+2][1]*F[1][k-i+2]: od: eqs := eqs union {temp}: od: for k from 2 to maxL+1 do temp := G[1][k]: for i from 1 to num do if Steps[i][1] > Steps[i][2] then temp := temp - z^Steps[i][2]*F[abs(Steps[i][1]-Steps[i][2])][k-1]: fi: od: eqs := eqs union {temp}: od: for k from 2 to maxL+1 do temp := G[k][1]: for i from 1 to num do if Steps[i][1] < Steps[i][2] then temp := temp - z^Steps[i][2]*F[k-1][abs(Steps[i][1]-Steps[i][2])]: fi: od: eqs := eqs union {temp}: od: vars := {}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do vars := vars union {F[i][j],G[i][j]}: od: od: vars := vars minus {F[1][1]}: sol := eliminate(eqs,vars): if nops({sol}) > 1 then return seq(subs(F[1][1]=F,sol[i][2]),i=1..nops({sol})): else return op(subs(F[1][1]=F,sol[2])): fi: end: #RGF2DWE(Steps,z,F,t) finds the rigorous generating function for an arbitrary sequence of steps with infinite width and variable t counting the number of times the walk touches the line y=x. RGF2DWE:= proc(Steps, z, H,t) local i,j,k,G,vars,eqs,maxL,temp,num,sol,diffop,recop: num := nops(Steps): maxL := 0: for i from 1 to num do if abs(Steps[i][1]-Steps[i][2]) > maxL+1 then maxL := abs(Steps[i][1]-Steps[i][2])-1: fi: od: eqs := F[1][1]-1-G[1][1]*F[1][1]: for i from 1 to nops(Steps) do if Steps[i][1] = Steps[i][2] then eqs := eqs-z^Steps[i][1]*F[1][1]: fi: od: eqs := {eqs}: temp := G[1][1]: for i from 1 to num do for j from 1 to num do if Steps[i][1] > Steps[i][2] and Steps[j][1] < Steps[j][2] then temp := temp - z^(Steps[j][2]+Steps[i][2])*F[abs(Steps[i][1]-Steps[i][2])][abs(Steps[j][1]-Steps[j][2])]: fi: od: od: eqs := eqs union {temp}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do if i <> 1 and j <> 1 then if i >= j then eqs := eqs union {G[i][j]-G[i-j+1][1]}: else eqs := eqs union {G[i][j]-G[1][j-i+1]}: fi: fi: od: od: for i from 1 to maxL+1 do for j from 1 to maxL+1 do if i <> j then if i <> 1 or j <> 1 then temp := F[i][j]: for k from 1 to min(i,j) do if i < j then temp := temp - F[i-k+1][1]*G[k][j]: else temp := temp - G[i][k]*F[1][j-k+1]: fi: od: eqs := eqs union {temp}: fi: fi: od: od: for k from 2 to maxL+1 do temp := F[k][k]-F[1][1]: for i from 2 to k do temp := temp - G[k-i+2][1]*F[1][k-i+2]: od: eqs := eqs union {temp}: od: for k from 2 to maxL+1 do temp := G[1][k]: for i from 1 to num do if Steps[i][1] > Steps[i][2] then temp := temp - z^Steps[i][2]*F[abs(Steps[i][1]-Steps[i][2])][k-1]: fi: od: eqs := eqs union {temp}: od: for k from 2 to maxL+1 do temp := G[k][1]: for i from 1 to num do if Steps[i][1] < Steps[i][2] then temp := temp - z^Steps[i][2]*F[k-1][abs(Steps[i][1]-Steps[i][2])]: fi: od: eqs := eqs union {temp}: od: temp := H-1: for i from 1 to nops(Steps) do if Steps[i][1] = Steps[i][2] then temp := temp-t*z^Steps[i][1]*H: fi: od: for i from 1 to num do for j from 1 to num do if Steps[i][1] > Steps[i][2] and Steps[j][1] < Steps[j][2] then temp := temp - t*z^(Steps[j][2]+Steps[i][2])*F[abs(Steps[i][1]-Steps[i][2])][abs(Steps[j][1]-Steps[j][2])]*H: fi: od: od: eqs := eqs union {temp}: vars := {}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do vars := vars union {F[i][j],G[i][j]}: od: od: sol := eliminate(eqs,vars): return op(sol[2]): if nops({sol}) > 1 then return seq(subs(sol[i],H),i=1..nops({sol})): else return subs(sol,H): fi: end: #rgfprove(Steps,ope,W,N) "proves" that ope(W,N) actually holds for arbitrary w by looking at the generating functions rgfprove:= proc(Steps,ope,W,N) local i,j,k,l,F,G,vars,eqs,maxL,temp,temp1,n,sol,sol2,ans,deg: #What the answer should look like: replay W^max ->1, W^(max-1)->F[1], W^(max-2) ->F[1]F[2], ... , etc. ans := 0: _EnvAllSolutions := true: for i from 1 to nops(ope) do deg := degree(op(i,ope),W): ans := ans + coeff(op(i,ope),W,deg)*mul(F[j][1][1],j=1..degree(ope,W)-deg): od: print(`The solution to be proved is:`,ans,`=0`): n := nops(Steps): maxL := 0: for i from 1 to n do if abs(Steps[i][1]-Steps[i][2]) > maxL+1 then maxL := abs(Steps[i][1]-Steps[i][2])-1: fi: od: eqs := F[1][1][1]-1-G[1][1][1]*F[1][1][1]: for i from 1 to nops(Steps) do if Steps[i][1] = Steps[i][2] then eqs := eqs-N^Steps[i][1]*F[1][1][1]: fi: od: eqs := {eqs}: temp := G[1][1][1]: for i from 1 to n do for j from 1 to n do if Steps[i][1] > Steps[i][2] and Steps[j][1] < Steps[j][2] then temp := temp - N^(Steps[j][2]+Steps[i][2])*F[2][abs(Steps[i][1]-Steps[i][2])][abs(Steps[j][1]-Steps[j][2])]: fi: od: od: eqs := eqs union {temp}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do if i <> 1 and j <> 1 then if i >= j then eqs := eqs union {G[1][i][j]-G[i][i-j+1][1]}: else eqs := eqs union {G[1][i][j]-G[j][1][j-i+1]}: fi: fi: od: od: for i from 1 to maxL+1 do for j from 1 to maxL+1 do if i <> j then if i <> 1 or j <> 1 then temp := F[1][i][j]: for k from 1 to min(i,j) do if i < j then temp := temp - F[1][i-k+1][1]*G[1][k][j]: else temp := temp - G[1][i][k]*F[1][1][j-k+1]: fi: od: eqs := eqs union {temp}: fi: fi: od: od: for k from 2 to maxL+1 do temp := F[1][k][k]-F[k][1][1]: for i from 2 to k do temp := temp - G[i-1][k-i+2][1]*F[i-1][1][k-i+2]: od: eqs := eqs union {temp}: od: for k from 2 to maxL+1 do temp := G[1][1][k]: for i from 1 to n do if Steps[i][1] > Steps[i][2] then if k-1 < abs(Steps[i][1]-Steps[i][2]) then temp := temp - N^Steps[i][2]*F[2][abs(Steps[i][1]-Steps[i][2])-k+2][1]: else temp := temp - N^Steps[i][2]*F[2][abs(Steps[i][1]-Steps[i][2])][k-1]: fi: fi: od: eqs := eqs union {temp}: od: for k from 2 to maxL+1 do temp := G[1][k][1]: for i from 1 to n do if Steps[i][1] < Steps[i][2] then if k-1 < abs(Steps[i][1]-Steps[i][2]) then temp := temp - N^Steps[i][2]*F[2][1][abs(Steps[i][1]-Steps[i][2])-k+2]: else temp := temp - N^Steps[i][2]*F[2][k-1][abs(Steps[i][1]-Steps[i][2])]: fi: fi: od: eqs := eqs union {temp}: od: vars := {}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do for k from 1 to maxL+1 do vars := vars union {G[k][i][j]}: od: od: od: sol := eliminate(eqs,vars)[2]: print(`The relation among various generating functions is`): print(sol): print(`with the following convention for F[i][j][k]: i stands for the symbolic width w-i+1`): print(`and j and k stand for the [j-1,k-1] walk`): sol2 := sol: for l from 1 to degree(ope,W)-1 do for i from 1 to nops(sol) do temp := 0: for j from 1 to nops(sol[i]) do temp1 := 1: if nops(op(j,sol[i])) > 1 then for k from 1 to nops(op(j,sol[i])) do if type(op(k,op(j,sol[i])),indexed) then temp1 := temp1*incby1(op(k,op(j,sol[i]))): else temp1 := temp1*op(k,op(j,sol[i])): fi: od: elif type(op(j,sol[i]),indexed) then temp1 := temp1*incby1(op(j,sol[i])): fi: temp := temp + temp1: od: sol2 := sol2 union {temp} od: sol := sol2 minus sol: od: vars := {}: for i from 1 to maxL+1 do for j from 1 to maxL+1 do for k from 1 to degree(ope,W) do vars := vars union {F[k][i][j]}: od: od: od: vars := vars minus {seq(F[i][1][1],i=1..degree(ope,W)),seq(seq(F[1][i][j],i=2..maxL+2),j=2..maxL+1)}: sol := solve(sol2): #return sol: print(`Using many copies of the above equations with different widths, solving for them`): print(`And substituting them in the expected solution gives:`): return {seq(simplify(subs(sol[i],ans)),i=1..nops({sol}))}; end: incby1 := proc(F) local i,j,k,temp: temp := op(0,op(0,op(0,F))): i := op(1,op(0,op(0,F)))+1: j := op(1,op(0,F)): k := op(1,F): return temp[i][j][k]: end: ##################Functions from GuessHolo2: findrec() and Yafe() #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,dim: option remember: dim := LinearAlgebra[Dimension](f); if (1+DEGREE)*(1+ORDER)+5+ORDER>dim then error "Insufficient data for a recurrence of order %1 degree %2",ORDER, DEGREE: fi: ope:=0: var:={}: ope := add( add(a[i,j]*n^j*N^i, j=0..DEGREE), i=0..ORDER); var := indets(ope) minus {n,N}; eq := {seq( add(subs(n=n0,coeff(ope,N,i))*f[n0+i], i=0..ORDER), n0 = 1..dim-ORDER)}; var1:=solve(eq,var): kv := map(lhs, select( x->evalb(op(1,x) = op(2,x)), var1)); ope:=eval(ope, var1); if ope=0 then return FAIL end if; ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: `if`( nops(ope)>=1,Yafe(ope[1],N)[2],FAIL); end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then return (1,0) end if; ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): # ope1:=normal(ope1): ope1 := collect(ope1, N, factor); factor(coe1),ope1: end: