###################################################################### ##DRUNKARD: Save this file as DRUNKARD. To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read DRUNKARD: # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Dec. 10, 2007 print(`Created: Dec. 10, 2007`): print(` This is DRUNKARD, `): print(`to enumerate the number of simple lattice walks`): print(`in (a) Z^k (b) x1>=0, ..., xk>=0, (c) x1>=x2>= ...>=xk `): print(` (d) x1>=x2>= ...>=xk>=0 `): print(`and to estimate the probability of returning to the origin`): print(`while staying at that region. `): print(`It accompanies the paper `): print(` How Likely Is Polya's Drunkard to Return From The Pub`): print(` Without Getting Mugged? (In d-Dimensional Manhattan)`): print(`by Doron Zeilberger`): print(`published in the Personal Journal of Ekhad and Zeilberger `): print(`http://www.math.rutgers.edu/~zeilberg/pj.html .`): print(``): print(`Please report bugs to: zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/drunkard.html .`): print(`For a list of the different helps, and for the MAIN procedures, type ezra();`): print(``): with(combinat): ezra:=proc(): if nargs=0 then print(`The MAIN procedures are: POLYA, POLYAc, SIPUR , `): print(`For help with it, type: ezra(SIPUR); ezra(POLYA), ezra(POLYAc)); `): print(`The help for the other procedures is arranged as follows. `): print(): print(`There are the following kinds of Help`): print(`ezra0: Unrestricted simple (and lazy) random walks in Z^k .`): print(`ezra1: Simple (and lazy) random walks in x1>=0, x2>=0, ..., xk>=0 .`): print(`ezra2: Simple (and lazy) random walks in x1>=x2>= ... >=xk .`): print(`ezra3: Simple (and lazy) random walks in x1>=x2>= ... >=xk>=0 .`): print(`ezraY: Help about counting Young Tableaux with <=k rows .`): print(`ezraGZ: Help about the Gessel-Zeilberger Reflection method .`): print(`ezraGF: Help about the Generating functions . `): print(`To get the list of relevant procedures in ezra0 do, ezra0(); `): print(`To get help of a procedure in ezra0, do ezra0(Procedurename)`): print(`Similarly for all the other ezras .`): elif nops([args])=1 and op(1,[args])=POLYA then print(`POLYA(K): all the Polya constants for walks`): print(`in (a) Z^k, (b) x1>=0, ..., xk>=0`): print(`(c) x1>=x2>= ...>=xk (d) x1>=x2>= ...>=xk>=0`): print(`for k=2..K`): print(`For example, try:`): print(`POLYA(3):`): elif nops([args])=1 and op(1,[args])=POLYAc then print(`POLYAc(K): Terse form of POLYA, given as a list of`): print(`quartets for k=2..K`): print(`all the Polya constants for walks`): print(`in (a) Z^k, (b) x1>=0, ..., xk>=0`): print(`(c) x1>=x2>= ...>=xk (d) x1>=x2>= ...>=xk>=0`): print(`for k=2..K`): print(`For example, try:`): print(`POLYAc(3):`): elif nops([args])=1 and op(1,[args])=SIPUR then print(`SIPUR(K,L,L1,n,N): all the stories for walks`): print(`in (a) Z^k, (b) x1>=0, ..., xk>=0`): print(`(c) x1>=x2>= ...>=xk (d) x1>=x2>= ...>=xk>=0`): print(`for k=1..K.`): print(`For example, try:`): print(`SIPUR(4,100,2000,n,N):`): else print(`There is no ezra for`, args): fi: end: ezraY:=proc() if args=NULL then print(`The procedures concerning Standard Young Tableaux are `): print(`Pars, Pkn, SeqPk, SeqYk, SeqYkr, SYT, Ykn, Yknr `): elif nops([args])=1 and op(1,[args])=DGB then print(`DGB(n): The Dominque Gouyou-Beauchamps formula for`): print(`the number of Young-Tableaux with n cells and <=5 rows.`): print(`For example, try: DGB(10);`): elif nops([args])=1 and op(1,[args])=Pars then print(`Pars(n,k): the set of partitions of n with <=k parts`): elif nops([args])=1 and op(1,[args])=SeqPk then print(`SeqPk(k,N): The sequence for the `): print(` number of permutations with of length n without an `): print(`increasing subsequence of length>k for n=1.. N. `): print(`For example, try:`): print(`SeqPk(3,10);`); elif nops([args])=1 and op(1,[args])=SeqYk then print(`SeqYk(k,N): The sequence for the `): print(` number of Standard Young Tableaux with n cells and`): print(`<=k rows, for n=1.. N. For example, try:`): print(`SeqYk(3,10);`); elif nops([args])=1 and op(1,[args])=SeqYkr then print(`SeqYkr(k,N,r): The sequence for sum of (f_lam )^r`): print(` over all partions of n with`): print(`<=k parts, for n=1.. N. For example, try:`): print(`SeqYkr(3,10,3);`); elif nops([args])=1 and op(1,[args])=SYT then print(`SYT(lam): the number of Standard Young Tableaux of shape lam.`): print(`For example, try: SYT([4,2,2]);`): elif nops([args])=1 and op(1,[args])=Ykn then print(`Ykn(k,n): The number of Standard Young Tableaux with n cells and`): print(`<=k rows. For example, try:`): print(`Ykn(3,10);`); elif nops([args])=1 and op(1,[args])=Pkn then print(`Pkn(k,n): The number of permutations of {1,..,n} `): print(` with no increasing subsequence of length>k. For example, try`): print(`Pkn(3,10);`); print(`There is no Young-Tableaux ezra (ezraY) for`,args): fi: end: ezraGZ:=proc() if args=NULL then print(`The procedures about the Gessel-Zeilbeger Reflection `): print(`Method are: `): print(` ASx0, ASx0x1, AS, Phi0, Phi1, Phi2, Phi3, Trans, Trans1 `): elif nops([args])=1 and op(1,[args])=AS then print(`AS(S,k,t,x): The antisymmetrizer of the monomial `): print(`t[1]^x[1]*...*t[k]^x[k] under the Coxeter group`): print(`generated by the set of reflections in S. For example, try:`): print(`AS({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x);`): elif nops([args])=1 and op(1,[args])=ASx0 then print(`ASx0(S,k,t,x,x0): The anti-symmetrizer of the monomial `): print(`t[1]^x0[1]*...*t[k]^x0[k] under the Coxeter group`): print(`generated by the set of reflections in S, where`): print(`x0 is a specific numerical point in Z^k.`): print(` For example, try:`): print(`ASx0({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x,[2,1,0]);`): elif nops([args])=1 and op(1,[args])=ASx0x1 then print(`ASx0x1(S,k,t,x,x0,x1): The anti-symmetrizer of the monomial `): print(`t[1]^x0[1]*...*t[k]^x0[k] under the Coxeter group`): print(`generated by the set of reflections in S, where`): print(`x0 is a specific numerical point in Z^k.`): print(`divided by t[1]^x1[1]* ... *t[k]^xk[k] `): print(` For example, try:`): print(`ASx0x1({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x,[2,1,0],[1,1,1]);`): elif nops([args])=1 and op(1,[args])=Phi0 then print(`Phi0(t,k): The Phi-polynomial for type-0 walks`): print(`which is just 1. Try: `): print(`Phi0(t,2);`): elif nops([args])=1 and op(1,[args])=Phi1 then print(`Phi1(t,k): The Phi-polynomial for type 1 walks`): print(`For example, try:`): print(`Phi1(t,2);`): elif nops([args])=1 and op(1,[args])=Phi2 then print(`Phi2(t,k): The Phi-polynomial for type 2 walks`): print(`For example, try:`): print(`Phi2(t,2);`): elif nops([args])=1 and op(1,[args])=Phi3 then print(`Phi3(t,k): The Phi-polynomial for type 3 walks`): print(`For example, try:`): print(`Phi3(t,2);`): elif nops([args])=1 and op(1,[args])=Trans then print(`Trans(S,k,x): Given a set S of affine-linear transformations`): print(`in x[1], ..., x[k] (x is a symbol and k is a positive integer)`): print(`and a non-negative integer d, finds all transformations generated`): print(`by S. For example, try:`): print(`Trans({[x[2],x[1]]} ,2,x);`): elif nops([args])=1 and op(1,[args])=Trans1 then print(`Trans1(S,k,x,d): Given a set S of affine-linear transformations`): print(`in x[1], ..., x[k] (x is a symbol and k is a positive integer)`): print(`and a non-negative integer d, finds all transformations whose`): print(`(minimal) length is d. For example, try:`): print(`Trans1({[x[2],x[1]]} ,2,x,1);`): else print(`There is no ezraGZ for`,args): fi: end: ezraGF:=proc() if args=NULL then print(`The procedures about the Generating Function Approach are: `): print(` JatN, MonoToJ, PolyToJ `): elif nops([args])=1 and op(1,[args])=JatN then print(`JatN(a,t,N): The first N terms of Sum(t^(2*n+a)/n!/(n+a)!, n=0..N);`): print(`For example, try: JatN(1,t,10);`): elif nops([args])=1 and op(1,[args])=MonoToJ then print(`MonoToJ(M,t,k,J): Given a monomial M in t[1], ..., t[k]`): print(`t[1]^a1*...*t[k]^ak, and a symbol J outputs J[a1]*...J[ak]`): print(`For example, try:`): print(`MonoToJ(t[1]^2*t[2],t,2,J);`); elif nops([args])=1 and op(1,[args])=PolyToJ then print(`PolyToJ(P,t,k,J): Given a polynomial P in t[1], ..., t[k]`): print(`and a symbol J outputs the umbral image t^a->J[a]`): print(`For example, try:`): print(`PolyToJ(t[1]+t[2],t,2,J);`): else print(`There is no ezraGZ for`,args): fi: end: ezra0:=proc() if args=NULL then print(`The procedures for Type 0 walks ( walks in Z^k ) `): print(`with no restrictions are `): print(` DrunkW0, ExV0, f0kn, g0kn, NuW0, GF0, GFb0, NuLW0,`): print(`Pr0, Seq0slow, Seq0, Seq0o`): print(` SeqD0, SeqD0t, SeqW0, SeqLW0 , `): print(` f0kn, g0kn, NuW0CT, NuLW0CT , PollW0, Polya0, Sipur0 `): print(): print(): print(` The MAIN procedures are: Seq0, GF0, Pr0, Sipur0 `): elif nops([args])=1 and op(1,[args])=D0kn then print(`D0kn(k,n): The number of`): print(`random walks starting and ending at the origin`): print(`in Z^k `): print(`without restrictions `): print(`taking n days, where at each day one has the option either to perform`): print(`no step, one step , or two steps in the unit directions.`): print(`For example, try: D0kn(2,10);`): elif nops([args])=1 and op(1,[args])=D0tkn then print(`D0tkn(k,n): The number of`): print(`random walks starting and ending at the origin`): print(`in Z^k `): print(`without restrictions `): print(`taking n days, where at each day one has the option either to perform`): print(`one step or two steps in the unit directions.`): print(`For example, try: D0tkn(2,10);`): elif nops([args])=1 and op(1,[args])=DrunkW0 then print(`DrunkW0(S,F,k,N): Walking from S to F in Z^k`): print(`true if it comes to F in <=N steps, false otherwise. `): print(`For example, try:`): print(`DrunkW0([0,0],[0,0],2,10000);`): elif nops([args])=1 and op(1,[args])=ExV0 then print(`ExV0(S,F): The expected number of visits of a random`): print(`walker from S to F in Z^k (k=nops(S)). For example`): print(` try: ExV0([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=f0kn then print(`f0kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n steps, in Z^k`): print(` using unit steps in`): print(`all directions. For example, try:`): print(`f0kn(2,10);`): elif nops([args])=1 and op(1,[args])=g0kn then print(`g0kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n days, in`): print(`Z^k, and where at at day you either walk one unit step`): print(`in some direction, or do nothing. For example, try:`): print(`g0kn(2,10);`): elif nops([args])=1 and op(1,[args])=NuW0 then print(`NuW0(S,F,n): The number of walks in Z^k with`): print(`unit steps in all directions, from pt. S to pt. F`): print(`having n steps. For example, try:`): print(`NuW0([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuW0CT then print(`NuW0CT(S,F,n): The number of walks in Z^k with`): print(`unit steps in all directions, from pt. S to pt. F`): print(`having n steps. Using the Constant Term formula. For example, try:`): print(`NuW0CT([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuLW0 then print(`NuLW0(S,F,n): The number of walks in Z^k with`): print(`unit steps in all directions, from pt. S to pt. F`): print(`as well as the lazy option of doing nothing, `): print(`having n steps. For example, try:`): print(`NuLW0([0,0],[5,3],11);`): elif nops([args])=1 and op(1,[args])=PollW0 then print(`PollW0(S,F,k,N,K): letting the walker walk from S to F in Z^k`): print(`with N=infinity, and doing it K times. Return the percentage`): print(`of returns. For example, try:`): print(`PollW0([0,0,0],[0,0,0],3,10000,100);`): elif nops([args])=1 and op(1,[args])=Polya0 then print(`Polya0(k,L,L1): The prob. of returning to the origin`): print(` in simple random walk in Z^k.`): print(`For example, try:`): print(`Polya0(3,50,1000);`): elif nops([args])=1 and op(1,[args])=Pr0 then print(`Pr0(S,F): The probability of a Polya drunkard staring`): print(`at S and making it to F walking in Z^k (k=nops(S)). `): print(`if S=F, that probabilty of visiting S at least one more time`): print(`For example try: Pr0([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=SeqD0 then print(`SeqD0(k,L): the sequence whose i^th term is D0kn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD0(2,20);`): elif nops([args])=1 and op(1,[args])=SeqD0t then print(`SeqD0t(k,L): the sequence whose i^th term is D0tkn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD0t(2,20);`): elif nops([args])=1 and op(1,[args])=SeqLW0 then print(`SeqLW0(k,L): the sequence whose i^th term is f0kn(k,i) `): print(` For example, try: `): print(`SeqLW0(2,20);`): elif nops([args])=1 and op(1,[args])=SeqW0 then print(`SeqW0(k,L): the sequence whose i^th term (for i=1..L), is the number`): print(`of simple random walks in Z^k, from the origin back`): print(`For example, try:`): print(`SeqW0(2,20);`): elif nops([args])=1 and op(1,[args])=Seq0slow then print(`Seq0slow(S,k,F,L): `): print(`the sequence whose i^th term (for i=1..L), is the number`): print(`of simple random walks in Z^k, S to F .`): print(`For example, try:`): print(`Seq0slow([0,0],[1,2], 2,20);`): elif nops([args])=1 and op(1,[args])=Sipur0 then print(`Sipur0(k,L,L1,n,N): The story of simple lattice walks in dimension`): print(`k, listing L`): print(`terms in the sequence. later using L1 terms to estimate the `): print(` Polya constant, i.e. the probability of return to the origin `): print(`For example, try:`): print(`Sipur0(2,50,1000, n,N);`): elif nops([args])=1 and op(1,[args])=GF0 then print(`GF0(S,F,J): The formal power series`): print(`(as a polynomial in J[a]'s, whose coeff. of t^n/n!`): print(`is the number of simple random walks`): print(`in Z^k `): print(`in F.It also outputs the largest index of J `): print(`For example, try:`): print(` GF0([0,0],[5,3],J); `): elif nops([args])=1 and op(1,[args])=GFb0 then print(`GFb0(S,F,x): GF0(S,F,x) in terms of the`): print(`Modified Bessel Functions BesselI `): print(`For example, try:`): print(` GFb0([0,0],[5,3],x); `): elif nops([args])=1 and op(1,[args])=Seq0 then print(`Seq0(S,F,L): Like Seq0slow but using (exp.) generating function`): print(`For example, try: Seq0([0,0],[0,0],20);`): elif nops([args])=1 and op(1,[args])=Seq0o then print(`Seq0o(k,L): The sequence, of length L, whose i-th`): print(`term is the number of walks, with unit steps in all directions`): print(`in Z^k`): print(`starting and ending at the origin `): print(`For example, try: Seq0o(2,10);`): else print(`There is no ezra0 for`,args): fi: end: ezra1:=proc() if args=NULL then print(`The procedures for Type 1 walks ( walks in Z^k ) `): print(`that stay in x1>=0, x2>=0, ..., xk>=0, are `): print(` f1kn, ExV1, g1kn, PollW1, Polya1`): print(` Pr1, NuW1, GF1, GFb1, NuW1CT, NuLW1, NuLW1CT,`): print(` Seq1slow, Seq1, Seq1o,SeqD1, SeqD1t, SeqW1, SeqLW1 , Polya1, Sipur1 `): print(): print(): print(` The MAIN procedures are: Seq1, GF1, Pr1, Sipur1 `): elif nops([args])=1 and op(1,[args])=D1kn then print(`D1kn(k,n): The number of`): print(`random walks starting and ending at the origin`): print(`in Z^k `): print(`staying in x1>=0, ..., xk>=0 `): print(`taking n days, where at each day one has the option either to perform`): print(`no step, one step , or two steps in the unit directions.`): print(`For example, try: D1kn(2,10);`): elif nops([args])=1 and op(1,[args])=D1tkn then print(`D1tkn(k,n): The number of`): print(`random walks starting and ending at the origin`): print(`in Z^k `): print(`staying in x1>=0, ..., xk>=0 `): print(`taking n days, where at each day one has the option either to perform`): print(`one step or two steps in the unit directions.`): print(`For example, try: D1tkn(2,10);`): elif nops([args])=1 and op(1,[args])=ExV1 then print(`ExV1(S,F): The expected number of visits of a random`): print(`walker from S to F in Z^k (k=nops(S)). `): print(`Staying in x1>=x2>= ...>=xk`): print(`For example`): print(`try: ExV1([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=f1kn then print(`f1kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n steps, in Z^k`): print(`staying in x1>=0, ..., xk>=0, `): print(` using unit steps in`): print(`all directions. For example, try:`): print(`f1kn(2,10);`): elif nops([args])=1 and op(1,[args])=g1kn then print(`g1kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n days, in`): print(`Z^k, `): print(`staying in x1>=0, ..., xk>=0 `): print(` and where at at day you either walk one unit step`): print(`in some direction, or do nothing. For example, try:`): print(`g1kn(2,10);`): elif nops([args])=1 and op(1,[args])=NuW1 then print(`NuW1(S,F,n): The number of walks in Z^k `): print(`staying in x1>=0, ..., xk>=0 `): print(`with unit steps in all directions, from pt. S to pt. F`): print(`having n steps. For example, try:`): print(`NuW1([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuW1CT then print(`NuW1CT(S,F,n): Like NuW1 but using the Constant Term formula`): print(`For example, try:`): print(`NuW1CT([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuLW1 then print(`NuLW1(S,F,n): The number of walks in Z^k `): print(`staying in x1>=0, ..., xk>=0 ,`): print(`with unit steps in all directions, from pt. S to pt. F`): print(`as well as the lazy option of doing nothing, `): print(`having n steps. For example, try:`): print(`NuLW1([0,0],[5,3],11);`): elif nops([args])=1 and op(1,[args])=PollW1 then print(`PollW1(S,F,k,N,K): letting the walker walk from S to F in Z^k`): print(`staying in x1>=0, ..., xk>=0 `): print(`with N=infinity, and doing it K times. Return the percentage`): print(`of returns. For example, try:`): print(`PollW1([0,0,0],[0,0,0],3,10000,100);`): elif nops([args])=1 and op(1,[args])=Polya1 then print(`Polya1(k,L,L1): The prob. of returning to the origin`): print(` in simple random walk in Z^k`): print(`staying in x1>=0, ..., xk>=0 `): print(`For example, try:`): print(`Polya1(3,50,1000);`): elif nops([args])=1 and op(1,[args])=Pr1 then print(`Pr1(S,F): The probability of a Polya drunkard staring`): print(`at S and making it to F walking in Z^k (k=nops(S)). `): print(`staying in x1>=0, ..., xk>=0 `): print(`if S=F, that probabilty of visiting S at least one more time`): print(`For example try: Pr1([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=SeqD1 then print(`SeqD1(k,L): the sequence whose i^th term is D1kn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD1(2,20);`): elif nops([args])=1 and op(1,[args])=SeqD1t then print(`SeqD1t(k,L): the sequence whose i^th term is D1tkn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD1t(2,20);`): elif nops([args])=1 and op(1,[args])=SeqLW1 then print(`SeqLW1(k,L): the sequence whose i^th term is f1kn(k,i) `): print(` For example, try: `): print(`SeqLW1(2,20);`): elif nops([args])=1 and op(1,[args])=SeqW1 then print(`SeqW1(k,L): the sequence whose i^th term (i=1..L), is the number`): print(`of simple random walks in Z^k, `): print(`staying in x1>=0, ..., xk>=0, `): print(` from the origin back to the origin .` ): print(`For example, try:`): print(`SeqW1(2,20);`): elif nops([args])=1 and op(1,[args])=Seq1slow then print(`Seq1slow(S,F,k,L): `): print(`the sequence whose i^th term (for i=1..L), is the number`): print(`of simple random walks `): print(`staying in x1>=0, ..., xk>=0, `): print(`from S to F`): print(`For example, try:`): print(`Seq1slow([0,0],[1,2], 2,20);`): elif nops([args])=1 and op(1,[args])=Sipur1 then print(`Sipur1(k,L,L1,n,N): The story of simple lattice walks in dimension`): print(`k, staying in x1>=0, ..., xk>=0 `): print(`listing L`): print(`terms in the sequence. later using L1 terms to estimate the `): print(` relative Polya constant. `): print(` the probability of return to the origin staying in the above region`): print(`For example, try:`): print(`Sipur1(2,50,1000, n,N);`): elif nops([args])=1 and op(1,[args])=GF1 then print(`GF1(S,F,J): The formal power series`): print(`(as a polynomial in J[a]'s, whose coeff. of t^n/n!`): print(`is the number of simple random walks`): print(`staying in x[1]>=0, ..., x[k]>=0 starting at S ending `): print(`in F.It also outputs the largest index of J `): print(`For example, try:`): print(` GF1([0,0],[5,3],J); `): elif nops([args])=1 and op(1,[args])=GFb1 then print(`GFb1(S,F,x): GF1(S,F,x) in terms of the`): print(`Modified Bessel Functions BesselI `): print(`For example, try:`): print(` GFb1([0,0],[5,3],x); `): elif nops([args])=1 and op(1,[args])=Seq1 then print(`Seq1(S,F,L): Like Seq1slow but using (exp.) generating function`): print(`For example, try: Seq1([0,0],[0,0],20);`): elif nops([args])=1 and op(1,[args])=Seq1o then print(`Seq1o(k,L): The sequence, of length L, whose i-th`): print(`term is the number of walks, with unit steps in all directions`): print(`in x1>0, x2>=0, ..., xk>=0 `): print(`starting and ending at the origin `): print(`For example, try: Seq1o(2,10);`): else print(`There is no ezra1 for`,args): fi: end: ezra2:=proc() if args=NULL then print(`The procedures for Type 1 walks ( walks in Z^k ) `): print(`that stay in x1>=x2>= ... >= xk , are `): print(` f2kn, ExV2`): print(` g2kn, NuW2, NuW2CT, GF2, GFb2, NuLW2, PollW2,`): print(` Polya2, Pr2, Seq2slow,Seq2 `): print(`Seq2o, SeqD2, SeqD2t, SeqW2, SeqL2 , Sipur2`): print(): print(): print(` The MAIN procedures are: Seq2, GF2, Pr2, Sipur2 `): elif nops([args])=1 and op(1,[args])=D2kn then print(`D2kn(k,n): The number of`): print(`random walks starting and ending at the origin`): print(`in Z^k `): print(`staying in x1>= x2>= ... >= xk `): print(`taking n days, where at each day one has the option either to perform`): print(`no step, one step , or two steps in the unit directions.`): print(`For example, try: D2kn(2,10);`): elif nops([args])=1 and op(1,[args])=D2tkn then print(`D2tkn(k,n): The number of`): print(`random walks starting and ending at the origin`): print(`in Z^k `): print(`staying in x1>= x2>= ... >= xk `): print(`taking n days, where at each day one has the option either to perform`): print(`one step or two steps in the unit directions.`): print(`For example, try: D2tkn(2,10);`): elif nops([args])=1 and op(1,[args])=ExV2 then print(`ExV2(S,F): The expected number of visits of a random`): print(`walker from S to F in Z^k (k=nops(S)). `): print(`Staying in x1>=x2>= ...>=xk`): print(`For example`): print(`try: ExV2([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=f2kn then print(`f2kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n steps, in Z^k`): print(`staying in x1>= x2>= ... >= xk `): print(` using unit steps in`): print(`all directions. For example, try:`): print(`f2kn(2,10);`): elif nops([args])=1 and op(1,[args])=g2kn then print(`g2kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n days, in`): print(`Z^k, `): print(`staying in x1>= x2>= ... >= xk `): print(` and where at at day you either walk one unit step`): print(`in some direction, or do nothing. For example, try:`): print(`g2kn(2,10);`): elif nops([args])=1 and op(1,[args])=NuW2 then print(`NuW2(S,F,n): The number of walks in Z^k `): print(`staying in x1>= x2>= ... >= xk `): print(`with unit steps in all directions, from pt. S to pt. F`): print(`having n steps. For example, try:`): print(`NuW2([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuW2CT then print(`NuW2CT(S,F,n): Like NuW2)S,F,n), but using the Constant Term `): print(`Formula . `): print(`For example, try:`): print(`NuW2CT([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuLW2 then print(`NuLW2(S,F,n): The number of walks in Z^k `): print(`staying in x1>= x2>= ... >= xk `): print(`with unit steps in all directions, from pt. S to pt. F`): print(`as well as the lazy option of doing nothing, `): print(`having n steps. For example, try:`): print(`NuLW2([0,0],[5,3],11);`): elif nops([args])=1 and op(1,[args])=PollW2 then print(`PollW2(S,F,k,N,K): letting the walker walk from S to F in Z^k`): print(`staying in x1>=x2>= ... >= xk `): print(`with N=infinity, and doing it K times. Return the percentage`): print(`of returns. For example, try:`): print(`PollW2([0,0,0],[0,0,0],3,10000,100);`): elif nops([args])=1 and op(1,[args])=Polya2 then print(`Polya2(k,L,L1): The prob. of returning to the origin`): print(` in simple random walk in Z^k`): print(`staying in x1>=x2>= ... xk `): print(`For example, try:`): print(`Polya2(3,50,1000);`): elif nops([args])=1 and op(1,[args])=Pr2 then print(`Pr2(S,F): The probability of a Polya drunkard staring`): print(`at S and making it to F walking in Z^k (k=nops(S)). `): print(`staying in x1>=x2>=...>=xk `): print(`if S=F, that probabilty of visiting S at least one more time`): print(`For example try: Pr2([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=SeqD2 then print(`SeqD2(k,L): the sequence whose i^th term is D2kn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD2(2,20);`): elif nops([args])=1 and op(1,[args])=SeqD2t then print(`SeqD2t(k,L): the sequence whose i^th term is D2tkn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD2t(2,20);`): elif nops([args])=1 and op(1,[args])=SeqLW2 then print(`SeqLW2(k,L): the sequence whose i^th term is f2kn(k,i) `): print(` For example, try: `): print(`SeqLW2(2,20);`): elif nops([args])=1 and op(1,[args])=SeqW2 then print(`SeqW2(k,L): the sequence whose i^th term (i=1..L), is the number`): print(`of simple random walks in Z^k, `): print(`staying in x1>=x2>= ...>= xk , `): print(` from the origin back`): print(`For example, try:`): print(`SeqW2(2,20);`): elif nops([args])=1 and op(1,[args])=Seq2slow then print(`Seq2slow(S,F,k,L): `): print(`the sequence whose i^th term (for i=1..L), is the number`): print(`of simple random walks `): print(`staying in x1>=x2>= ...>= xk , `): print(`from S to F`): print(`For example, try:`): print(`Seq2slow([0,0],[1,2], 2,20);`): elif nops([args])=1 and op(1,[args])=Sipur2 then print(`Sipur2(k,L,L1,n,N): The story of simple lattice walks in dimension`): print(`k, staying in x1>=x2> ... >= xk `): print(`listing L`): print(`terms in the sequence. later using L1 terms to estimate the `): print(` relative Polya constant, i.e. `): print(` the probability of return to the origin staying in the above region`): print(`For example, try:`): print(`Sipur2(2,50,1000, n,N);`): elif nops([args])=1 and op(1,[args])=GF2 then print(`GF2(S,F,J): The formal power series`): print(`(as a polynomial in J[a]'s, whose coeff. of t^n/n!`): print(`is the number of simple random walks`): print(`staying in x[1]>=...>=x[k] starting at S ending `): print(`in F.It also outputs the largest index of J `): print(`For example, try:`): print(` GF2([0,0],[5,3],J); `): elif nops([args])=1 and op(1,[args])=GFb2 then print(`GFb2(S,F,x): GF2(S,F,x) in terms of the`): print(`Modified Bessel Functions BesselI `): print(`For example, try:`): print(` GFb2([0,0],[5,3],x); `): elif nops([args])=1 and op(1,[args])=Seq2 then print(`Seq2(S,F,L): Like Seq2slow but using (exp.) generating function`): print(`For example, try: Seq2([0,0],[0,0],20);`): elif nops([args])=1 and op(1,[args])=Seq2o then print(`Seq2o(k,L): The sequence, of length L, whose i-th`): print(`term is the number of walks, with unit steps in all directions`): print(`in x1>=x2>= ...>=xk `): print(`starting and ending at the origin `): print(`For example, try: Seq2o(2,10);`): else print(`There is no ezra2 for`,args): fi: end: ezra3:=proc() if args=NULL then print(`The procedures for Type 3 walks (walks that stay in:`): print(`x1>=x2>= ...>xk>=0, are: `): print(` f3kn, ExV3, g3kn, NuW3, NuW3CT, GF3, GFb3, `): print(` NuLW3, PollW3, Polya3, Pr3, Seq3slow, Seq3, Seq3o, `): print(`SeqD3, SeqD3t, SeqW3, SeqLW3 , `): print(): print(): print(` The MAIN procedures are: Seq3, GF3, Pr3, Sipur3 `): elif nops([args])=1 and op(1,[args])=D3kn then print(`D3kn(k,n): The number of (k+1)-noncrossing tangled diagrams over [n]`): print(`also the number of random walks starting at and ending at the origin`): print(`staying in the region x1>=x2>= ...>=xk>=0 `): print(`taking n days, where at each day one has the option either to perform`): print(`no step, one step , or two steps in the unit directions.`): print(`For example, try: D3kn(2,10);`): elif nops([args])=1 and op(1,[args])=D3tkn then print(`D3tkn(k,n): The number of (k+1)-noncrossing tangled diagrams over [n]`): print(`also the number of random walks starting at and ending at the origin`): print(`taking n days, where at each day one has the option either to perform`): print(`one or two random steps.`): print(`For example, try: D3tkn(2,10);`): elif nops([args])=1 and op(1,[args])=ExV3 then print(`ExV3(S,F): The expected number of visits of a random`): print(`walker from S to F in Z^k (k=nops(S)). `): print(`Staying in x1>=x2>= ...>=xk>=0`): print(`For example`): print(`try: ExV3([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=f3kn then print(`f3kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n steps, staying`): print(`in x1>=x2>=...>=xk>=0 and using unit steps in`): print(`all directions. For example, try:`): print(`f3kn(2,10);`): elif nops([args])=1 and op(1,[args])=g3kn then print(`g3kn(k,n): The number of ways of going from the`): print(`origin back to the origin in n days, staying`): print(`in x1>=x2>=...>=xk>=0 and where at at day you walk one unit step`): print(`in some direction, or do nothing. For example, try:`): print(`g3kn(2,10);`): elif nops([args])=1 and op(1,[args])=NuW3 then print(`NuW3(S,F,n): The number of walks in x1>=x2>=...>=xk>=0 with`): print(`unit steps in all directions, from pt. S to pt. F`): print(`having n steps. For example, try:`): print(`NuW3([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuW3CT then print(`NuW3CT(S,F,n): Like NuW3(S,F,n), but using the constant term`): print(`formula `): print(`For example, try:`): print(`NuW3CT([0,0],[5,3],12);`): elif nops([args])=1 and op(1,[args])=NuLW3 then print(`NuLW3(S,F,n): The number of walks in x1>=x2>=...>=xk>=0 with`): print(`unit steps in all directions, from pt. S to pt. F`): print(`as well as the lazy option of doing nothing, `): print(`having n steps. For example, try:`): print(`NuLW3([0,0],[5,3],11);`): elif nops([args])=1 and op(1,[args])=GF3 then print(`GF3(S,F,J): The formal power series`): print(`(as a polynomial in J[a]'s, whose coeff. of t^n/n!`): print(`is the number of simple random walks`): print(`staying in x[1]>=...>=x[k]>=0 starting at S ending `): print(`in F.It also outputs the largest index of J `): print(`For example, try:`): print(` GF3([0,0],[5,3],J); `): elif nops([args])=1 and op(1,[args])=GFb3 then print(`GFb3(S,F,x): GF3(S,F,x) in terms of the`): print(`Modified Bessel Functions BesselI `): print(`For example, try:`): print(` GFb2([0,0],[5,3],x); `): elif nops([args])=1 and op(1,[args])=PollW3 then print(`PollW3(S,F,k,N,K): letting the walker walk from S to F in Z^k`): print(`staying in x1>=x2>= ... >= xk >=0 `): print(`with N=infinity, and doing it K times. Return the percentage`): print(`of returns. For example, try:`): print(`PollW3([0,0,0],[0,0,0],3,10000,100);`): elif nops([args])=1 and op(1,[args])=Polya3 then print(`Polya3(k,L,L1): The prob. of returning to the origin`): print(` in simple random walk in Z^k`): print(`staying in x1>=x2>= ... xk>=0 `): print(`For example, try:`): print(`Polya3(3,50,1000);`): elif nops([args])=1 and op(1,[args])=Pr3 then print(`Pr3(S,F): The probability of a Polya drunkard staring`): print(`at S and making it to F walking in Z^k (k=nops(S)). `): print(`staying in x1>=x2>=...>=xk `): print(`if S=F, that probabilty of visiting S at least one more time`): print(`For example try: Pr3([0,0,0],[0,0,0]);`): elif nops([args])=1 and op(1,[args])=Seq3 then print(`Seq3(S,F,L): Like Seq3 but using (exp.) generating function`): print(`For example, try: Seq3([0,0],[0,0],20);`): elif nops([args])=1 and op(1,[args])=Seq3o then print(`Seq3o(k,L): The sequence, of length L, whose i-th`): print(`term is the number of walks, with unit steps in all directions`): print(`in x1>=x2>= ...>=xk>=0 `): print(`starting and ending at the origin `): print(`For example, try: Seq3o(2,10);`): elif nops([args])=1 and op(1,[args])=SeqD3 then print(`SeqD3(k,L): the sequence whose i^th term is D3kn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD3(2,20);`): elif nops([args])=1 and op(1,[args])=SeqD3t then print(`SeqD3t(k,L): the sequence whose i^th term is D3tkn(k,i) `): print(`For i=1..L. For example, try:`): print(`SeqD3t(2,20);`): elif nops([args])=1 and op(1,[args])=SeqLW3 then print(`SeqLW3(k,L): the sequence whose i^th term is the number`): print(`of lazy simple random walks from the origin back`): print(` to the origin with`): print(`i steps staying at x1>=x2>=...>=xk>=0. For example, try:`): print(`SeqLW3(2,20);`): elif nops([args])=1 and op(1,[args])=SeqW3 then print(`SeqW3(k,L): the sequence whose i^th term is the number`): print(`of simple random walks from the origin back`): print(` to the origin with`): print(`i steps staying at x1>=x2>=...>=xk>=0. For example, try:`): print(`SeqW3(2,20);`): elif nops([args])=1 and op(1,[args])=Seq3slow then print(`Seq3slow(S,F,k,L): `): print(`the sequence whose i^th term (for i=1..L), is the number`): print(`of simple random walks `): print(`staying in x1>=x2>= ...>=xk<=0 , `): print(`from S to F`): print(`For example, try:`): print(`Seq3slow([0,0],[1,2], 2,20);`): elif nops([args])=1 and op(1,[args])=Sipur3 then print(`Sipur3(k,L,L1,n,N): The story of simple lattice walks in dimension`): print(`k, staying in x1>=x2> ... >= xk>=0 `): print(`listing L`): print(`terms in the sequence. later using L1 terms to estimate the `): print(` relative Polya constant (the prob. of returning to the origin `): print(`and staying in the above region `): print(`For example, try:`): print(`Sipur3(2,50,1000, n,N);`): else print(`There is no ezra3 for`,args): fi: end: ######Begin SYT##### #DGB(n): The Dominque Gouyou-Beauchamps formula for #the number of Young-Tableaux with n cells and <=5 rows. #For example, try: DGB(10); DGB:=proc(n) local i: add(3!*n!*(2*i+2)!/(n-2*i)!/i!/(i+1)!/(i+2)!/(i+3)!,i=0..trunc(n/2)): end: #SYT(lam): the number of Standard Young Tableaux of Shape lam #For example, try: SYT([4,2,2]); SYT:=proc(lam) local i,gu,k: option remember: if lam=[] then RETURN(1): fi: k:=nops(lam): gu:=0: for i from 1 to k-1 do if lam[i]>lam[i+1] then gu:=gu+SYT([op(1..i-1,lam),lam[i]-1,op(i+1..k,lam)]): fi: od: if lam[k]>1 then gu:=gu+SYT([op(1..k-1,lam),lam[k]-1]): else gu:=gu+SYT([op(1..k-1,lam)]): fi: gu: end: #Pars(n,k): the set of partitions of n with <=k parts Pars:=proc(n,k) local k1: {seq(op(Pars1(n,k1)),k1=1..k)}: end: #Pars1(n,k): the set of integer-partitions of n with exactly k parts Pars1:=proc(n,k) local gu,g,i1: option remember: if nk then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: else add(NuW0(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuW0(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k): fi: end: #NuW0CT(S,F,n): The number of walks in Z^k with #unit steps in all directions, from pt. S to pt. F #having n steps. #Using the Constant-Term formula #For example, try: #NuW0CT([0,0],[5,3],5); NuW0CT:=proc(S,F,n) local k,i,gu,x,lu: k:=nops(S): gu:=add(x[i]+1/x[i],i=1..k): gu:=expand(gu^n): lu:=mul(x[i]^S[i],i=1..k): gu:=expand(gu*lu): for i from 1 to k do gu:=coeff(gu,x[i],F[i]): od: gu: end: #CheckNuW0CT(S,N,T): checks that NuW0 and NuW0CT coincide #for starting point S and all points in Pars(N,k) and time t<=T #For example, try: CheckNuW0CT([3,2],7,10); CheckNuW0CT:=proc(S,N,T) local k,gu,t,j: k:=nops(S): gu:=Pars(N,k) minus Pars(N,k-1): {seq(seq(evalb(NuW0(S,gu[j],t)=NuW0CT(S,gu[j],t)), j=1..nops(gu)), t=1..T)}: end: #f0kn(k,n): The number of ways of going from the #origin back to the origin in n steps, in Z^k, with no restrictions #using unit steps in #all directions. For example, try: #f0kn(2,10); f0kn:=proc(k,n): NuW0([0$k],[0$k],n) :end: #g0kn(k,n): The number of ways of going from the #origin back to the origin in n steps, in Z^k #unit steps in #all directions plus the lazy option. For example, try: #g0kn(2,10); g0kn:=proc(k,n) local i: add(binomial(n,i)*f0kn(k,n-i),i=0..n) :end: #SeqW0(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #2*i steps For example, try: #SeqW0(2,20); SeqW0:=proc(k,L) local i: [seq(f0kn(k,2*i),i=1..L)]:end: Seq0slow:=proc(S,F,k,L) local i: [seq(NuW0(S,F,i),i=1..L)]:end: #GF0(S,F,J): The formal power series #(as a polynomial in J[a]'s, whose coeff. of t^n/n! #is the number of simple random walks #in F. It also outputs the largest index of J #For example, try: #GF0([0,0],[5,3],J); GF0:=proc(S,F,J) local k,i,t,lu,gadol: k:=nops(S): lu:=Phi0G(t,S): lu:=expand(lu/mul(t[i]^F[i],i=1..k)): gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)): PolyToJ(lu,t,k,J),gadol: end: GFb0:=proc(S,F,x) local gu,J,L,i: gu:=GF0(S,F,J): L:=gu[2]: gu:=gu[1]: subs({seq(J[i]=BesselI(i,2*x),i=0..L)}, gu): end: #ExV0(S,F): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)). For example #try: ExV0([0,0,0],[0,0,0]); ExV0:=proc(S,F) local x,gu,d: option remember: d:=nops(S): if nops(F)<>d then ERROR(`bad input`): fi: gu:=GFb0(S,F,x): evalf(Int(subs(x=x/(2*d),gu)*exp(-x),x=0..infinity)): end: #ExV0d(S,F,m): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)),doing it directly, #walking up to m steps #. For example #try: ExV0d([0,0,0],[0,0,0],100); ExV0d:=proc(S,F,m) local x,gu,d,k: k:=nops(S): gu:=Seq0(S,F,m): if S=F then 1.+add(gu[i]/(2.*k)^i,i=1..nops(gu)): else add(gu[i]/(2.*k)^i,i=1..nops(gu)): fi: end: #Pr0(S,F): The probability of a Polya drunkard staring #at S and making it to F walking in Z^k (k=nops(S)). #if S=F, that probabilty of visiting S at least one more time #For example try: Pr0([0,0,0],[0,0,0]); Pr0:=proc(S,F) local q,gu: if S=F then gu:=ExV0(S,F): RETURN(1-1/gu): fi: q:=Pr0(F,F): ExV0(S,F)*(1-q): end: #Seq0(S,F,L): Like Seq0slow but using (exp.) generating function #For example, try: Seq0([0,0],[0,0],2,20); Seq0:=proc(S,F,L) local k,gu,J,gadol,t,i: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: gu:=GF0(S,F,J): gadol:=gu[2]: gu:=gu[1]: for i from 0 to gadol do gu:=subs(J[i]=JatN(i,t,L),gu): od: gu:=taylor(gu,t=0,L+1): [seq(coeff(gu,t,i)*i!,i=1..L)]: end: Seq0o:=proc(k,L) local gu,i: gu:=Seq0([0$k],[0$k],2*L): [seq(gu[2*i],i=1..L)]: end: #SeqW0gf(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #2*i steps. Using the generating function expression: For example, try: #SeqW0gf(2,20); SeqW0gf:=proc(k,L) local i,gu,t: gu:=add(t^(2*i)/i!^2,i=0..L): gu:=gu^k: gu:=taylor(gu,t=0,2*L+1): [seq(coeff(gu,t,2*i)*(2*i)!,i=1..L)]: end: #NuLW0(S,F,n): The number of walks in Z^k #unit steps in all directions, from pt. S to pt. F #as well as the lazy (do-nothing) step #having n steps. For example, try: #NuLW0([0,0],[5,3],5); NuLW0:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: else add(NuLW0(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuLW0(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k)+ NuLW0(S,F,n-1): : fi: end: #SeqLW0Old(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. For example, try: #SeqLW0Old(2,20); SeqLW0Old:=proc(k,L) local i: [seq(NuLW0([0$k],[0$k],i),i=1..L)]:end: #SeqLW0(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. For example, try: #SeqLW0(2,20); SeqLW0:=proc(k,L) local i: [seq(g0kn(k,i),i=1..L)]:end: #SeqLW0gf(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. Using the generating function. For example, try: #SeqLW0gf(2,20); SeqLW0gf:=proc(k,L) local gu,i,n,t: gu:=add(t^(2*i)/i!^2,i=0..L): gu:=gu^k: gu:=taylor(gu,t=0,2*L+1): [seq(add(binomial(n,i)*(n-i)!*coeff(gu,t,n-i),i=0..n),n=1..L)]: end: #D0tkn(k,n): The number of walks taking n days in Z^k #where at each day, you can either walk one steps or two steps #For example, try: D0tkn(2,10); D0tkn:=proc(k,n) local i: add(binomial(n,i)*f0kn(k,2*n-i),i=0..n):end: #SeqD0t(k,L): The sequence of D0tkn(k,i) for i=1..L. For example, try: #SeqD0t(2,20); SeqD0t:=proc(k,L) local i: [seq(D0tkn(k,i),i=1..L)]:end: #D0kn(k,n):The number of walks taking n days in Z^k #where at each day, you can either walk zero steps, one steps or two steps #For example, try: D0kn(2,10); D0kn:=proc(k,n) local i: add(binomial(n,i)*D0tkn(k,n-i),i=0..n):end: #SeqD0(k,L): The sequence of D0kn(k,i) for i=1..L. For example, try: #SeqD0(2,20); SeqD0:=proc(k,L) local i: [seq(D0kn(k,i),i=1..L)]:end: #Sipur0(k,L,L1,n,N): The story of simple lattice walks in dimension #k, and listing L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Sipur0(2,50,1000,n,N); Sipur0:=proc(k,L,L1,n,N) local gu,C,ope,lu,Kavua,J,a,t: print(`This is the story of simple random walks that start and`): print(`end at the origin, in the`, k, `-dimensional (cubic) lattice`): print(`The exponential Generating Function is`): print(GF0([0$k],[0$k],J)[1]): print(`where J[a](t) is: `): print(Sum(t^(2*n+a)/n!/(n+a)!,n=0..infinity)): gu:=Seq0o(k,L): print(`The sequence of the number of such walks of length 2n for`): print(` n=1..`, L, `is: `): print(gu): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then print(L , `terms do not suffice to guess a recurrence `): fi: print(`The sequence is annihilated by the recurrence operator`): print(ope): if ope<>FAIL then gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): fi: if ope<>FAIL then lu:=Asy(ope,n,N,3): print(`The asymptotic behaviour is`): Kavua:=evalf( gu[nops(gu)]/subs(n=nops(gu),lu) ): print(Kavua*lu) : fi: if k>=3 then print(`The expected number of visits to the origin is `): print(ExV0([0$k],[0$k])): print(`The Polya constant, i.e. the probablity of returning to the origin`): print(`is estimated to be`): print(Pr0([0$k],[0$k])): fi: end: #Polya0(k,L,L1): The Polya constant in dim k #For example, try: #Polya0(2,50,1000,n,N); Polya0:=proc(k,L,L1) local gu,C,ope,m,err,Di,exp1,n,N,i: if k<=2 then RETURN(1): fi: gu:=Seq0o(k,L): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then print(L , `terms do not suffice to guess a recurrence `): fi: if ope<>FAIL then gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): fi: m:=1.0+add(evalf(gu[i]/(2*k)^(2*i)),i=1..nops(gu)): exp1:=1-k/2: err:=(nops(gu))^exp1: Di:=min(trunc(-ln(err)/ln(10))-3,10); if Di>2 then evalf(1-1/m,Di): fi: end: #DrunkW0(S,F,k,N): Walking from S to F in Z^k #true if it comes to F in <=N steps, false otherwise #For example, try: #DrunkW0([0,0],[0,0],2,10000); DrunkW0:=proc(S,F,k,N) local ra,x,i,r: if {nops(S),nops(F)}<>{k} then RETUNR(FAIL): fi: x:=S: ra:=rand(1..2*k): for i from 1 to N do r:=ra(): if r mod 2 = 0 then r:=r/2: x:=[op(1..r-1,x),x[r]+1,op(r+1..k,x)]: else r:=(r+1)/2: x:=[op(1..r-1,x),x[r]-1,op(r+1..k,x)]: fi: if x=F then RETURN(true): fi: od: false: end: #PollW0(S,F,k,N,K): letting the walker walk from S to F in Z^k #with N=infinity, and doing it K times. Return the percentage #of returns. For example, try: #PollW0([0,0,0],[0,0,0],3,10000,100); PollW0:=proc(S,F,k,N,K) local c,i: c:=0: for i from 1 to K do if DrunkW0(S,F,k,N) then c:=c+1: fi: od: evalf(c/K): end: #####End Walks in Z^k without restrictions #####Begin Walks in Z^k that stay in x1>=0, ..., xk>=0 #NuW1(S,F,n): The number of walks in Z^k with #unit steps in all directions, #staying in x1>=0, ..., xk>=0, from pt. S to pt. F, #having n steps. For example, try: #NuW1([0,0],[5,3],5); NuW1:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: elif {seq(evalb(F[i]>=0),i=1..k)}<>{true} then 0: else add(NuW1(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuW1(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k): fi: end: #NuW1CT(S,F,n): The number of walks in Z^k with #unit steps in all directions, staying in x[1]>=0 ... x[k]>=0 #from pt. S to pt. F #having n steps. #Using the Constant-Term formula #For example, try: #NuW1CT([0,0],[5,3],5); NuW1CT:=proc(S,F,n) local k,i,gu,t,lu: k:=nops(S): gu:=add(t[i]+1/t[i],i=1..k): gu:=expand(gu^n): lu:=Phi1G(t,S): gu:=expand(gu*lu): for i from 1 to k do gu:=coeff(gu,t[i],F[i]): od: gu: end: #CheckNuW1CT(S,N,T): checks that NuW1 and NuW1CT coincide #for starting point S and all points in Pars(N,k) and time t<=T #For example, try: CheckNuW1CT([3,2],7,10); CheckNuW1CT:=proc(S,N,T) local k,gu,t,j: k:=nops(S): gu:=Pars(N,k) minus Pars(N,k-1): {seq(seq(evalb(NuW1(S,gu[j],t)=NuW1CT(S,gu[j],t)), j=1..nops(gu)), t=1..T)}: end: #f1kn(k,n): The number of ways of going from the #origin back to the origin in n steps, in Z^k, #staying in x1>=0, ..., xk>=0 #using unit steps in #all directions. For example, try: #f1kn(2,10); f1kn:=proc(k,n): NuW1([0$k],[0$k],n) :end: #g1kn(k,n): The number of ways of going from the #origin back to the origin in n steps, in Z^k #staying in x1>=0, ..., xk>=0 #unit steps in #all directions plus the lazy option. For example, try: #g1kn(2,10); g1kn:=proc(k,n) local i: add(binomial(n,i)*f1kn(k,n-i),i=0..n) :end: #SeqW1(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #2*i steps, staying in x1>=0, ..., xk>=0 . For example, try: #SeqW1(2,20); SeqW1:=proc(k,L) local i: [seq(f1kn(k,2*i),i=1..L)]:end: Seq1slow:=proc(S,F,k,L) local i: [seq(NuW1(S,F,i),i=1..L)]:end: #GF1(S,F,J): The formal power series #(as a polynomial in J[a]'s, whose coeff. of t^n/n! #is the number of simple random walks #staying in x[1]>=0, ..., x[k]>=0 starting at S ending #in F. It also outputs the largest index of J #For example, try: #GF1([0,0],[5,3],J); GF1:=proc(S,F,J) local k,i,t,lu,gadol: k:=nops(S): lu:=Phi1G(t,S): lu:=expand(lu/mul(t[i]^F[i],i=1..k)): gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)): PolyToJ(lu,t,k,J),gadol: end: GFb1:=proc(S,F,x) local gu,J,L,i: gu:=GF1(S,F,J): L:=gu[2]: gu:=gu[1]: subs({seq(J[i]=BesselI(i,2*x),i=0..L)}, gu): end: #ExV1(S,F): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)). #Staying in x1>=0, ..., xk>=0 #For example #try: ExV1([0,0,0],[0,0,0]); ExV1:=proc(S,F) local x,gu,d: option remember: d:=nops(S): if nops(F)<>d then ERROR(`bad input`): fi: gu:=GFb1(S,F,x): evalf(Int(subs(x=x/(2*d),gu)*exp(-x),x=0..infinity)): end: #ExV1d(S,F,m): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)),region 1, doing it directly, #walking up to m steps #. For example #try: ExV1d([0,0,0],[0,0,0],100); ExV1d:=proc(S,F,m) local x,gu,d,k: k:=nops(S): gu:=Seq1(S,F,m): if S=F then 1.+add(gu[i]/(2.*k)^i,i=1..nops(gu)): else add(gu[i]/(2.*k)^i,i=1..nops(gu)): fi: end: #Pr1(S,F): The probability of a Polya drunkard staring #at S and making it to F walking in Z^k (k=nops(S)). #Staying in x1>=0, ..., xk>=0 #if S=F, that probabilty of visiting S at least one more time #For example try: Pr1([0,0,0],[0,0,0]); Pr1:=proc(S,F) local q,gu: if S=F then gu:=ExV1(S,F): RETURN(1-1/gu): fi: q:=Pr1(F,F): ExV1(S,F)*(1-q): end: #Seq1(S,F,L): Like Seq1slow but using (exp.) generating function #For example, try: Seq1([0,0],[0,0],20); Seq1:=proc(S,F,L) local k,gu,J,gadol,t,i: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: gu:=GF1(S,F,J): gadol:=gu[2]: gu:=gu[1]: for i from 0 to gadol do gu:=subs(J[i]=JatN(i,t,L),gu): od: gu:=taylor(gu,t=0,L+1): [seq(coeff(gu,t,i)*i!,i=1..L)]: end: Seq1o:=proc(k,L) local gu,i: gu:=Seq1([0$k],[0$k],2*L): [seq(gu[2*i],i=1..L)]: end: #NuLW1(S,F,n): The number of walks in Z^k, #staying in x1>=0, ..., xk>=0 #unit steps in all directions, from pt. S to pt. F #as well as the lazy (do-nothing) step #having n steps. For example, try: #NuLW1([0,0],[5,3],5); NuLW1:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: elif {seq(evalb(F[i]>=0),i=1..k)}<>{true} then 0: else add(NuLW1(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuLW1(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k)+ NuLW1(S,F,n-1): : fi: end: #SeqLW1Old(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. #staying in x1>=0, ..., xk>=0 #For example, try: #SeqLW1Old(2,20); SeqLW1Old:=proc(k,L) local i: [seq(NuLW1([0$k],[0$k],i),i=1..L)]:end: #SeqLW1(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. #staying in x1>=0, ..., xk>=0. For example, try: #SeqLW1(2,20); SeqLW1:=proc(k,L) local i: [seq(g1kn(k,i),i=1..L)]:end: #D1tkn(k,n): The number of walks taking n days in Z^k #staying in x1>=0, ..., xk>=0 #where at each day, you can either walk one steps or two steps #For example, try: D1tkn(2,10); D1tkn:=proc(k,n) local i: add(binomial(n,i)*f1kn(k,2*n-i),i=0..n):end: #SeqD1t(k,L): The sequence of D1tkn(k,i) for i=1..L. For example, try: #SeqD1t(2,20); SeqD1t:=proc(k,L) local i: [seq(D1tkn(k,i),i=1..L)]:end: #D1kn(k,n):The number of walks taking n days in Z^k #staying in x1>=0, ..., xk>=0 #where at each day, you can either walk zero steps, one steps or two steps #For example, try: D1kn(2,10); D1kn:=proc(k,n) local i: add(binomial(n,i)*D1tkn(k,n-i),i=0..n):end: #SeqD1(k,L): The sequence of D1kn(k,i) for i=1..L. For example, try: #SeqD1(2,20); SeqD1:=proc(k,L) local i: [seq(D1kn(k,i),i=1..L)]:end: #Sipur1(k,L,L1,n,N): The story of simple lattice walks in dimension #k, staying in x1>=0, ..., xk>=0, and listing L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Sipur1(2,50,1000,n,N); Sipur1:=proc(k,L,L1,n,N) local gu,C,ope,lu,Kavua,J,a,t,i,x: print(`This is the story of simple random walks that start and`): print(`end at the origin, in the`, k, `-dimensional (cubic) lattice`): print([seq(x[i]>=0, i=1..k)]): print(`The exponential Generating Function is`): print(GF1([0$k],[0$k],J)[1]): print(`where J[a](t) is: `): print(Sum(t^(2*n+a)/n!/(n+a)!,n=0..infinity)): gu:=Seq1o(k,L): print(`The sequence of the number of such walks of length 2n for`): print(` n=1..`, L, `is: `): print(gu): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then print(L , `terms do not suffice to guess a recurrence `): fi: print(`The sequence is annihilated by the recurrence operator`): print(ope): if ope<>FAIL then gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): fi: if ope<>FAIL then lu:=Asy(ope,n,N,3): print(`The asymptotic behaviour is`): Kavua:=evalf( gu[nops(gu)]/subs(n=nops(gu),lu) ): print(Kavua*lu) : fi: print(`The expected number of visits to the origin is `, ExV1([0$k],[0$k])): print(`The relative Polya constant, i.e. the probablity of `): print(`returning to the origin, while staying in `): print([seq(x[i]>=0, i=1..k)]): print(`is `, Pr1([0$k],[0$k])): end: #Polya1(k,L,L1): The Polya const. for #k, staying in x1>=0, ..., xk>=0, and listing L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Polya1(2,50,1000); Polya1:=proc(k,L,L1) local gu,C,ope,m,err,exp1,Di,n,N,i: gu:=Seq1o(k,L): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then RETURN(FAIL): fi: gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): exp1:=1-3*k/2: err:=(nops(gu))^exp1: Di:=min(trunc(-ln(err)/ln(10))-3,10); if Di>=2 then m:=1.0+add(evalf(gu[i]/(2*k)^(2*i)),i=1..nops(gu)): evalf(1-1/m,Di): fi: end: #DrunkW1(S,F,k,N): Walking from S to F in Z^k #in x1>=0, ..., xk>=0 #true if it comes to F in <=N steps, false otherwise #For example, try: #DrunkW0([0,0],[0,0],2,10000); DrunkW1:=proc(S,F,k,N) local ra,x,i,r: if {nops(S),nops(F)}<>{k} then RETUNR(FAIL): fi: x:=S: ra:=rand(1..2*k): for i from 1 to N do r:=ra(): if r mod 2 = 0 then r:=r/2: x:=[op(1..r-1,x),x[r]+1,op(r+1..k,x)]: else r:=(r+1)/2: x:=[op(1..r-1,x),x[r]-1,op(r+1..k,x)]: fi: if min(op(x))<0 then RETURN(false): fi: if x=F then RETURN(true): fi: od: false: end: #PollW1(S,F,k,N,K): letting the walker walk from S to F in Z^k #staying in x1>=0, ..., xk>=0 #with N=infinity, and doing it K times. Return the percentage #of returns. For example, try: #PollW1([0,0,0],[0,0,0],3,10000,100); PollW1:=proc(S,F,k,N,K) local c,i: c:=0: for i from 1 to K do if DrunkW1(S,F,k,N) then c:=c+1: fi: od: evalf(c/K): end: #####End Walks in Z^k that stay in x1>=0, ..., xk>=0 #####Begin Walks in Z^k that stay in x1>=x2>= ...>=xk #NuW2(S,F,n): The number of walks in Z^k with #unit steps in all directions, #staying in x1>=x2>=...>=xk #having n steps. For example, try: #NuW2([0,0],[5,3],5); NuW2:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: elif {seq(evalb(F[i]>=F[i+1]),i=1..k-1)}<>{true} then 0: else add(NuW2(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuW2(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k): fi: end: #NuW2CT(S,F,n): The number of walks in Z^k with #unit steps in all directions, staying in x[1]>=x[2]>= ...>=x[k] #from pt. S to pt. F #having n steps. #Using the Constant-Term formula #For example, try: #NuW2CT([0,0],[5,3],5); NuW2CT:=proc(S,F,n) local k,i,gu,t,lu: k:=nops(S): gu:=add(t[i]+1/t[i],i=1..k): gu:=expand(gu^n): lu:=Phi2G(t,S): gu:=expand(gu*lu): for i from 1 to k do gu:=coeff(gu,t[i],F[i]): od: gu: end: #CheckNuW2CT(S,N,T): checks that NuW2 and NuW2CT coincide #for starting point S and all points in Pars(N,k) and time t<=T #For example, try: CheckNuW2CT([3,2],7,10); CheckNuW2CT:=proc(S,N,T) local k,gu,t,j: k:=nops(S): gu:=Pars(N,k) minus Pars(N,k-1): [seq(seq(NuW2(S,gu[j],t)-NuW2CT(S,gu[j],t), j=1..nops(gu)), t=1..T)]: end: #f2kn(k,n): The number of ways of going from the #origin back to the origin in n steps, in Z^k, #staying in x1>=x2>=...>=xk #using unit steps in #all directions. For example, try: #f2kn(2,10); f2kn:=proc(k,n): NuW2([0$k],[0$k],n) :end: #g2kn(k,n): The number of ways of going from the #origin back to the origin in n steps, in Z^k #staying in x1>=x2>=...>=xk #unit steps in #all directions plus the lazy option. For example, try: #g2kn(2,10); g2kn:=proc(k,n) local i: add(binomial(n,i)*f2kn(k,n-i),i=0..n) :end: #SeqW2(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #2*i steps, staying in #x1>=x2>= ... >=xk . For example, try: #SeqW2(2,20); SeqW2:=proc(k,L) local i: [seq(f2kn(k,2*i),i=1..L)]:end: Seq2slow:=proc(S,F,k,L) local i: [seq(NuW2(S,F,i),i=1..L)]:end: #GF2(S,F,J): The formal power series #(as a polynomial in J[a]'s, whose coeff. of t^n/n! #is the number of simple random walks #staying in x[1]>=...>=x[k] starting at S ending #in F. It also outputs the largest index of J #For example, try: #GF2([0,0],[5,3],J); GF2:=proc(S,F,J) local k,i,t,lu,gadol: k:=nops(S): lu:=Phi2G(t,S): lu:=expand(lu/mul(t[i]^F[i],i=1..k)): gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)): PolyToJ(lu,t,k,J),gadol: end: GFb2:=proc(S,F,x) local gu,J,L,i: gu:=GF2(S,F,J): L:=gu[2]: gu:=gu[1]: subs({seq(J[i]=BesselI(i,2*x),i=0..L)}, gu): end: #ExV2(S,F): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)). #Staying in x1>=x2>= ...>= xk #For example #try: ExV2([0,0,0],[0,0,0]); ExV2:=proc(S,F) local x,gu,d: option remember: d:=nops(S): if nops(F)<>d then ERROR(`bad input`): fi: gu:=GFb2(S,F,x): evalf(Int(subs(x=x/(2*d),gu)*exp(-x),x=0..infinity)): end: #ExV2d(S,F,m): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)),region 2, doing it directly, #walking up to m steps #. For example #try: ExV2d([0,0,0],[0,0,0],100); ExV2d:=proc(S,F,m) local x,gu,d,k: k:=nops(S): gu:=Seq2(S,F,m): if S=F then 1.+add(gu[i]/(2.*k)^i,i=1..nops(gu)): else add(gu[i]/(2.*k)^i,i=1..nops(gu)): fi: end: #Pr2(S,F): The probability of a Polya drunkard staring #at S and making it to F walking in Z^k (k=nops(S)). #Staying in x1>=x2>=...>= xk #if S=F, that probabilty of visiting S at least one more time #For example try: Pr2([0,0,0],[0,0,0]); Pr2:=proc(S,F) local q,gu: if S=F then gu:=ExV2(S,F): RETURN(1-1/gu): fi: q:=Pr2(F,F): ExV2(S,F)*(1-q): end: #Seq2(S,F,L): Like Seq2slow but using (exp.) generating function #For example, try: Seq2([0,0],[0,0],20); Seq2:=proc(S,F,L) local k,gu,J,gadol,t,i: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: gu:=GF2(S,F,J): gadol:=gu[2]: gu:=gu[1]: for i from 0 to gadol do gu:=subs(J[i]=JatN(i,t,L),gu): od: gu:=taylor(gu,t=0,L+1): [seq(coeff(gu,t,i)*i!,i=1..L)]: end: Seq2o:=proc(k,L) local gu,i: gu:=Seq2([0$k],[0$k],2*L): [seq(gu[2*i],i=1..L)]: end: #NuLW2(S,F,n): The number of walks in Z^k, #staying in x1>=x2>= ... >=xk . #unit steps in all directions, from pt. S to pt. F #as well as the lazy (do-nothing) step #having n steps. For example, try: #NuLW2([0,0],[5,3],5); NuLW2:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: elif {seq(evalb(F[i]>=F[i+1]),i=1..k-1)}<>{true} then 0: else add(NuLW2(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuLW2(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k)+ NuLW2(S,F,n-1): : fi: end: #SeqLW2Old(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. #staying in x1>=x2>= ... >=xk . #For example, try: #SeqLW2Old(2,20); SeqLW2Old:=proc(k,L) local i: [seq(NuLW2([0$k],[0$k],i),i=1..L)]:end: #SeqLW2(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps in Z^k. #staying in x1>=x2>= ... >=xk . #For example, try: #SeqLW2(2,20); SeqLW2:=proc(k,L) local i: [seq(g2kn(k,i),i=1..L)]:end: #D2tkn(k,n): The number of walks taking n days in Z^k #staying in x1>=x2>= ... >=xk . #where at each day, you can either walk one steps or two steps #For example, try: D1tkn(2,10); D2tkn:=proc(k,n) local i: add(binomial(n,i)*f2kn(k,2*n-i),i=0..n):end: #SeqD2t(k,L): The sequence of D2tkn(k,i) for i=1..L. For example, try: #SeqD2t(2,20); SeqD2t:=proc(k,L) local i: [seq(D2tkn(k,i),i=1..L)]:end: #D2kn(k,n):The number of walks taking n days in Z^ #staying in x1>=x2>= ... >=xk . #where at each day, you can either walk zero steps, one steps or two steps #For example, try: D2kn(2,10); D2kn:=proc(k,n) local i: add(binomial(n,i)*D2tkn(k,n-i),i=0..n):end: #SeqD2(k,L): The sequence of D2kn(k,i) for i=1..L. For example, try: #SeqD2(2,20); SeqD2:=proc(k,L) local i: [seq(D2kn(k,i),i=1..L)]:end: #Sipur2(k,L,L1,n,N): The story of simple lattice walks in dimension #k, staying in x1>=x2>= ... >=xk , and listing L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Sipur2(2,50,1000,n,N); Sipur2:=proc(k,L,L1,n,N) local gu,C,ope,lu,Kavua,J,a,t,i,x: print(`This is the story of simple random walks that start and`): print(`end at the origin, in the`, k, `-dimensional (cubic) lattice`): print([seq(x[i]>=x[i+1], i=1..k-1)]): print(`The exponential Generating Function is`): print(GF2([0$k],[0$k],J)[1]): print(`where J[a](t) is: `): print(Sum(t^(2*n+a)/n!/(n+a)!,n=0..infinity)): gu:=Seq2o(k,L): print(`The sequence of the number of such walks of length 2n for`): print(` n=1..`, L, `is: `): print(gu): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then print(L , `terms do not suffice to guess a recurrence `): fi: print(`The sequence is annihilated by the recurrence operator`): print(ope): if ope<>FAIL then gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): fi: if ope<>FAIL then lu:=Asy(ope,n,N,3): print(`The asymptotic behaviour is`): Kavua:=evalf( gu[nops(gu)]/subs(n=nops(gu),lu) ): print(Kavua*lu) : fi: if k>=2 then print(`The expected number of visits to the origin is `, ExV2([0$k],[0$k])): print(`The relative Polya constant, i.e. the probablity of `): print(`returning to the origin, while staying in `): print([seq(x[i]>=x[i+1], i=1..k-1)]): print(` is equal to `, Pr2([0$k],[0$k])): fi: end: #Polya2(k,L,L1): The relative Polya constant for #simple lattice walks in dimension #k, staying in x1>=x2>= ... >=xk , and using L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Polya2(2,50,1000); Polya2:=proc(k,L,L1) local gu,C,ope,m,err,exp1,Di,n,N,i: gu:=Seq2o(k,L): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then RETURN(FAIL): fi: gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): exp1:=1-k^2/2: err:=(nops(gu))^exp1: Di:=min(trunc(-ln(err)/ln(10))-2,10); if k>=2 and Di>=2 then m:=1.0+add(evalf(gu[i]/(2*k)^(2*i)),i=1..nops(gu)): evalf(1-1/m,Di): else FAIL: fi: end: #DrunkW2(S,F,k,N): Walking from S to F in Z^k #in x1>=x2> ... >=xk #true if it comes to F in <=N steps, false otherwise #For example, try: #DrunkW2([0,0],[0,0],2,10000); DrunkW2:=proc(S,F,k,N) local ra,x,i,r,i1: if {nops(S),nops(F)}<>{k} then RETUNR(FAIL): fi: x:=S: ra:=rand(1..2*k): for i from 1 to N do r:=ra(): if r mod 2 = 0 then r:=r/2: x:=[op(1..r-1,x),x[r]+1,op(r+1..k,x)]: else r:=(r+1)/2: x:=[op(1..r-1,x),x[r]-1,op(r+1..k,x)]: fi: if min( seq(x[i1]-x[i1+1],i1=1..nops(x)-1) )<0 then RETURN(false): fi: if x=F then RETURN(true): fi: od: false: end: #PollW2(S,F,k,N,K): letting the walker walk from S to F in Z^k #staying in x1>=x2>= ..., >=xk #with N=infinity, and doing it K times. Return the percentage #of returns. For example, try: #PollW2([0,0,0],[0,0,0],3,10000,100); PollW2:=proc(S,F,k,N,K) local c,i: c:=0: for i from 1 to K do if DrunkW2(S,F,k,N) then c:=c+1: fi: od: evalf(c/K): end: #####End Walks in Z^k that stay in x1>=x2>= ...>=xk #####Begin Walks in x1>=x2>=...xk>=0 #NuW3(S,F,n): The number of walks in x1>=x2>=...>=xk>=0 with #unit steps in all directions, from pt. S to pt. F #having n steps. For example, try: #NuW3([0,0],[5,3],5); NuW3:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: elif {seq(evalb(F[i]>=F[i+1]),i=1..k-1),evalb(F[k]>=0)}<>{true} then 0: else add(NuW3(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuW3(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k): fi: end: #NuW3CT(S,F,n): The number of walks in Z^k with #unit steps in all directions, staying in type 3 wedge #from pt. S to pt. F #having n steps. #Using the Constant-Term formula #For example, try: #NuW3CT([0,0],[5,3],5); NuW3CT:=proc(S,F,n) local k,i,gu,t,lu: k:=nops(S): gu:=add(t[i]+1/t[i],i=1..k): gu:=expand(gu^n): lu:=Phi3G(t,S): gu:=expand(gu*lu): for i from 1 to k do gu:=coeff(gu,t[i],F[i]): od: gu: end: #GF3(S,F,J): The formal power series #(as a polynomial in J[a]'s, whose coeff. of t^n/n! #is the number of simple random walks #staying in x[1]>=...>=x[k]>=0 starting at S ending #in F. It also outputs the largest index of J #For example, try: #GF3([0,0],[5,3],J); GF3:=proc(S,F,J) local k,i,gu,t,lu,gadol: k:=nops(S): lu:=Phi3G(t,S): lu:=expand(lu/mul(t[i]^F[i],i=1..k)): gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)): PolyToJ(lu,t,k,J),gadol: end: GFb3:=proc(S,F,x) local gu,J,L,i: gu:=GF3(S,F,J): L:=gu[2]: gu:=gu[1]: subs({seq(J[i]=BesselI(i,2*x),i=0..L)}, gu): end: #ExV3(S,F): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)). #Staying in x1>=x2>= ...>=xk>=0 #For example #try: ExV3([0,0,0],[0,0,0]); ExV3:=proc(S,F) local x,gu,d: option remember: d:=nops(S): if nops(F)<>d then ERROR(`bad input`): fi: gu:=GFb3(S,F,x): evalf(Int(subs(x=x/(2*d),gu)*exp(-x),x=0..infinity)): end: #ExV3d(S,F,m): The expected number of visits of a random #walker from S to F in Z^k (k=nops(S)),region 3, doing it directly, #walking up to m steps #. For example #try: ExV3d([0,0,0],[0,0,0],100); ExV3d:=proc(S,F,m) local x,gu,d,k: k:=nops(S): gu:=Seq3(S,F,m): if S=F then 1.+add(gu[i]/(2.*k)^i,i=1..nops(gu)): else add(gu[i]/(2.*k)^i,i=1..nops(gu)): fi: end: #Pr3(S,F): The probability of a Polya drunkard staring #at S and making it to F walking in Z^k (k=nops(S)). #Staying in x1>=x2>=...>= xk>=0 #if S=F, that probabilty of visiting S at least one more time #For example try: Pr3([0,0,0],[0,0,0]); Pr3:=proc(S,F) local q,gu: if S=F then gu:=ExV3(S,F): RETURN(1-1/gu): fi: q:=Pr3(F,F): ExV3(S,F)*(1-q): end: #CheckNuW3CT(S,N,T): checks that NuW3 and NuW3CT coincide #for starting point S and all points in Pars(N,k) and time t<=T #For example, try: CheckNuW3CT([3,2],7,10); CheckNuW3CT:=proc(S,N,T) local k,gu,t,j: k:=nops(S): gu:=Pars(N,k) minus Pars(N,k-1): [seq(seq(NuW3(S,gu[j],t)-NuW3CT(S,gu[j],t), j=1..nops(gu)), t=1..T)]: end: #GF3(S,F,J): The formal power series #(as a polynomial in J[a]'s, whose coeff. of t^n/n! #is the number of simple random walks #staying in x[1]>=...>=x[k]>=0 starting at S ending #in F. It also outputs the largest index of J #For example, try: #GF3([0,0],[5,3],J); GF3:=proc(S,F,J) local k,i,t,lu,gadol: k:=nops(S): lu:=Phi3G(t,S): lu:=expand(lu/mul(t[i]^F[i],i=1..k)): gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)): PolyToJ(lu,t,k,J),gadol: end: #Seq3(S,F,L): Like Seq3slow but using (exp.) generating function #For example, try: Seq3([0,0],[0,0],20); Seq3:=proc(S,F,L) local k,gu,J,gadol,t,i: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: gu:=GF3(S,F,J): gadol:=gu[2]: gu:=gu[1]: for i from 0 to gadol do gu:=subs(J[i]=JatN(i,t,L),gu): od: gu:=taylor(gu,t=0,L+1): [seq(coeff(gu,t,i)*i!,i=1..L)]: end: Seq3o:=proc(k,L) local gu,i: gu:=Seq3([0$k],[0$k],2*L): [seq(gu[2*i],i=1..L)]: end: #f3kn(k,n): The number of ways of going from the #origin back to the origin in n steps, staying #in x1>=x2>=...>=xk>=0 and using unit steps in #all directions. For example, try: #f3kn(2,10); f3kn:=proc(k,n): NuW3([0$k],[0$k],n) :end: #g3kn(k,n): The number of ways of going from the #origin back to the origin in n steps, staying #in x1>=x2>=...>=xk>=0 and using unit steps in #all directions plus the lazy option. For example, try: #g3kn(2,10); g3kn:=proc(k,n) local i: add(binomial(n,i)*f3kn(k,n-i),i=0..n) :end: #SeqW3(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #2*i steps staying at x1>=x2>=...>=xk>=0. For example, try: #SeqW3(2,20); SeqW3:=proc(k,L) local i: [seq(f3kn(k,2*i),i=1..L)]:end: #Seq3slow(S,F,k,L): #the sequence whose i^th term (for i=1..L), is the number #of simple random walks #staying in x1>=x2>= ...>=xk>=0 #from S to F`): #For example, try: #Seq3slow([0,0],[1,2], 2,20);`): Seq3slow:=proc(S,F,k,L) local i: [seq(NuW3(S,F,i),i=1..L)]:end: #NuLW3(S,F,n): The number of walks in x1>=x2>=...>=xk>=0 with #unit steps in all directions, from pt. S to pt. F #as well as the lazy (do-nothing) step #having n steps. For example, try: #NuLW3([0,0],[5,3],5); NuLW3:=proc(S,F,n) local k,i: option remember: k:=nops(S): if nops(F)<>k then RETURN(FAIL): fi: if n=0 then if F=S then 1: else 0: fi: elif {seq(evalb(F[i]>=F[i+1]),i=1..k-1),evalb(F[k]>=0)}<>{true} then 0: else add(NuLW3(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+ add(NuLW3(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k)+ NuLW3(S,F,n-1): : fi: end: #SeqLW3Old(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps staying at x1>=x2>=...>=xk>=0. For example, try: #SeqLW3Old(2,20); SeqLW3Old:=proc(k,L) local i: [seq(NuLW3([0$k],[0$k],i),i=1..L)]:end: #SeqLW3(k,L): the sequence whose i^th term is the number #of non-lazy simple random walks from the origin back to the origin with #i steps staying at x1>=x2>=...>=xk>=0. For example, try: #SeqLW3(2,20); SeqLW3:=proc(k,L) local i: [seq(g3kn(k,i),i=1..L)]:end: #D3tkn(k,n): The number of (k+1)-noncrossing tangled diagrams over [n] #For example, try: D3tkn(2,10); D3tkn:=proc(k,n) local i: add(binomial(n,i)*f3kn(k,2*n-i),i=0..n):end: #SeqD3t(k,L): The sequence of D3tkn(k,i) for i=1..L. For example, try: #SeqD3t(2,20); SeqD3t:=proc(k,L) local i: [seq(D3tkn(k,i),i=1..L)]:end: #D3kn(k,n): The number of (k+1)-tangled diagrams over [n] #For example, try: D3kn(2,10); D3kn:=proc(k,n) local i: add(binomial(n,i)*D3tkn(k,n-i),i=0..n):end: #SeqD3(k,L): The sequence of D3kn(k,i) for i=1..L. For example, try: #SeqD3(2,20); SeqD3:=proc(k,L) local i: [seq(D3kn(k,i),i=1..L)]:end: #Sipur3(k,L,L1,n,N): The story of simple lattice walks in dimension #k, staying in x1>=x2>= ... >=xk>=0 , and listing L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Sipur3(2,50,1000,n,N); Sipur3:=proc(k,L,L1,n,N) local gu,C,ope,lu,Kavua,J,a,t,i,x: print(`This is the story of simple random walks that start and`): print(`end at the origin, in the`, k, `-dimensional (cubic) lattice`): print([seq(x[i]>=x[i+1], i=1..k-1),x[k]>=0]): print(`The exponential Generating Function is`): print(GF3([0$k],[0$k],J)[1]): print(`where J[a](t) is: `): print(Sum(t^(2*n+a)/n!/(n+a)!,n=0..infinity)): gu:=Seq3o(k,L): print(`The sequence of the number of such walks of length 2n for`): print(` n=1..`, L, `is: `): print(gu): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then print(L , `terms do not suffice to guess a recurrence `): fi: print(`The sequence is annihilated by the recurrence operator`): print(ope): if ope<>FAIL then gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): fi: if ope<>FAIL then lu:=Asy(ope,n,N,3): print(`The asymptotic behaviour is`): Kavua:=evalf( gu[nops(gu)]/subs(n=nops(gu),lu) ): print(Kavua*lu) : fi: print(`The expected number of visits to the origin is `, ExV3([0$k],[0$k])): print(`The relative Polya constant, i.e. the probablity of `): print(`returning to the origin, while staying in `): print([seq(x[i]>=x[i+1], i=1..k-1),x[k]>=0]): print(`is equal to`, Pr3([0$k],[0$k])) : end: #Polya3(k,L,L1): The Polya constant in Z^k #staying in x1>=x2>= ... >=xk>=0 , and using L #terms in the sequence, later using L1 terms #to estimate the Polya constant. #For example, try: #Polya3(2,50,1000); Polya3:=proc(k,L,L1) local gu,C,ope,m,err,exp1,Di,n,N,i: gu:=Seq3o(k,L): C:=trunc(2*sqrt(L)-5): ope:=Findrec(gu,n,N,C): if ope=FAIL then RETURN(FAIL): fi: gu:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1): exp1:=1-k^2-k/2: err:=(nops(gu))^exp1: Di:=min(trunc(-ln(err)/ln(10))-3,10); if Di>=3 then m:=1.0+add(evalf(gu[i]/(2*k)^(2*i)),i=1..nops(gu)): evalf(1-1/m,Di): else FAIL: fi: end: #DrunkW3(S,F,k,N): Walking from S to F in Z^k #in x1>=x2> ... >=xk>=0 #true if it comes to F in <=N steps, false otherwise #For example, try: #DrunkW3([0,0],[0,0],2,10000); DrunkW3:=proc(S,F,k,N) local ra,x,i,r,i1: if {nops(S),nops(F)}<>{k} then RETUNR(FAIL): fi: x:=S: ra:=rand(1..2*k): for i from 1 to N do r:=ra(): if r mod 2 = 0 then r:=r/2: x:=[op(1..r-1,x),x[r]+1,op(r+1..k,x)]: else r:=(r+1)/2: x:=[op(1..r-1,x),x[r]-1,op(r+1..k,x)]: fi: if min( seq(x[i1]-x[i1+1],i1=1..nops(x)-1), x[k] )<0 then RETURN(false): fi: if x=F then RETURN(true): fi: od: false: end: #PollW3(S,F,k,N,K): letting the walker walk from S to F in Z^k #staying in x1>=x2>= ..., >=xk>=0 #with N=infinity, and doing it K times. Return the percentage #of returns. For example, try: #PollW3([0,0,0],[0,0,0],3,10000,100); PollW3:=proc(S,F,k,N,K) local c,i: c:=0: for i from 1 to K do if DrunkW3(S,F,k,N) then c:=c+1: fi: od: evalf(c/K): end: #####End Walks in x1>=x2>=...xk>=0 ###Begin Gessel-Zeilberger Formula#### #Trans1(S,k,x,d): Given a set S of affine-linear transformations #in x[1], ..., x[k] (x is a symbol and k is a positive integer) #and a non-negative integer d, finds all transformations whose #(minimal) length is d. For example, try: #Trans1({[x[2],x[1]]} ,2,x,1); Trans1:=proc(S,k,x,d) local gu,i,mu,lu,s,kha,m: option remember: if d=0 then RETURN({[seq(x[i],i=1..k)]}): fi: mu:=Trans1(S,k,x,d-1): lu:={seq(op(Trans1(S,k,x,i)),i=0..d-1)}: gu:={}: for m in mu do for s in S do kha:=subs({seq(x[i]=s[i],i=1..k)},m): if not member(kha,lu) then gu:=gu union {kha}: fi: od: od: gu: end: #Trans(S,k,x): Given a set S of affine-linear transformations #in x[1], ..., x[k] (x is a symbol and k is a positive integer) #and a non-negative integer d, finds all transformations generated #by S. For example, try: #Trans({[x[2],x[1]]} ,2,x); Trans:=proc(S,k,x) local gu,d: gu:={}: for d from 0 while Trans1(S,k,x,d)<>{} do gu:=gu union Trans1(S,k,x,d): od: gu: end: #AS(S,k,t,x): The anti-symmetrizer of the monomial #t[1]^x[1]*...*t[k]^x[k] under the Coxeter group #generated by the set of reflections in S. For example, try: #AS({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x); AS:=proc(S,k,t,x) local gu,d,i,m,mu: gu:=0: for d from 0 while Trans1(S,k,x,d)<>{} do mu:=Trans1(S,k,x,d): gu:=gu+(-1)^d*add(mul(t[i]^m[i],i=1..k), m in mu): od: gu: end: #ASx0(S,k,t,x,x0): The anti-symmetrizer of the monomial #t[1]^x0[1]*...*t[k]^x0[k] under the Coxeter group #generated by the set of reflections in S, where #x0 is a specific numerical point in Z^k. # For example, try: #ASx0({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x,[2,1,0]); ASx0:=proc(S,k,t,x,x0) local gu,i: gu:=AS(S,k,t,x): factor(subs({seq(x[i]=x0[i],i=1..k)},gu)): end: #ASx0x1(S,k,t,x,x0,x1): The anti-symmetrizer of the monomial #t[1]^x0[1]*...*t[k]^x0[k] under the Coxeter group #generated by the set of reflections in S, where #x0 is a specific numerical point in Z^k. #Divided by t[1]^x1[1]*...*t[k]^x1[k] # For example, try: #ASx0x1({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x,[2,1,0],[1,1,1]); ASx0x1:=proc(S,k,t,x,x0,x1) local gu,i: gu:=ASx0(S,k,t,x,x0): gu:=gu/mul(t[i]^x1[i],i=1..k): gu:=expand(gu): end: Phi0:=proc(t,k) : 1: end: #Phi0G(t,S): The monomial #For example, try: #Phi0G(t,[0,0]); Phi0G:=proc(t,S) local S1,x,i,j,k: k:=nops(S): S1:={}: ASx0(S1,k,t,x,S): end: #Phi1(t,k): The Phi-polynomial for the group of reflections #w.r.t. x[i]=-1 . For example, try: #Phi1(t,2); Phi1:=proc(t,k) local S,x,i,j: S:={seq([seq(x[j], j=1..i-1),-2-x[i],seq(x[j],j=i+1..k)],i=1..k)}: ASx0(S,k,t,x,[0$k]): end: #Phi1G(t,S): The Anti-symmetrixer of t^S for the group of reflections #w.r.t. x[i]=-1 . For example, try: #Phi1G(t,[0,0]); Phi1G:=proc(t,S) local S1,x,i,j,k: k:=nops(S): S1:={seq([seq(x[j], j=1..i-1),-2-x[i],seq(x[j],j=i+1..k)],i=1..k)}: ASx0(S1,k,t,x,S): end: #Phi1F(t,k): The Phi-polynomial for the group of reflections #w.r.t. x[i]=-1, according to the product formula . #For example, try: #Phi1F(t,2); Phi1F:=proc(t,k) local i: mul(1-1/t[i]^2,i=1..k):end: #Phi2(t,k): The Phi-polynomial for the group of reflections #w.r.t. x[i]-x[i+1]=-1 . For example, try: #Phi2(t,2); Phi2:=proc(t,k) local S,x,i,j: S:={seq([seq(x[j], j=1..i-1),x[i+1]-1,x[i]+1, seq(x[j],j=i+2..k)],i=1..k-1)}: ASx0(S,k,t,x,[0$k]): end: #Phi2G(t,S): The anti-symmetrizer of t^S for the group of reflections #w.r.t. x[i]-x[i+1]=-1 . For example, try: #Phi2G(t,[0,0]); Phi2G:=proc(t,S) local S1,x,i,j,k: k:=nops(S): S1:={seq([seq(x[j], j=1..i-1),x[i+1]-1,x[i]+1, seq(x[j],j=i+2..k)],i=1..k-1)}: ASx0(S1,k,t,x,S): end: #Phi2F(t,k): The Phi-polynomial for the group of reflections #w.r.t. x[i]-x[i+1]=-1, computed by the product formula. For example, try: Phi2F:=proc(t,k) local i,j: mul(mul(t[i]-t[j],j=i+1..k),i=1..k)/mul(t[i]^(k-i),i=1..k): end: #Phi3(t,k): The Phi-polynomial for the group of reflections #w.r.t. x[i]-x[i+1]=-1 and x[k]=-1 . For example, try: #Phi3(t,2); Phi3:=proc(t,k) local S,x,i,j: S:={seq([seq(x[j], j=1..i-1),x[i+1]-1,x[i]+1, seq(x[j],j=i+2..k)],i=1..k-1), [seq(x[j],j=1..k-1),-2-x[k]]}: ASx0(S,k,t,x,[0$k]): end: #Phi3G(t,S): The anti-symmetrizerof t^S for the group of reflections #w.r.t. x[i]-x[i+1]=-1 and x[k]=-1 . For example, try: #Phi3G(t,[0,0]); Phi3G:=proc(t,S) local S1,x,i,j,k: k:=nops(S): S1:={seq([seq(x[j], j=1..i-1),x[i+1]-1,x[i]+1, seq(x[j],j=i+2..k)],i=1..k-1), [seq(x[j],j=1..k-1),-2-x[k]]}: ASx0(S1,k,t,x,S): end: #Phi3F(t,k): The Phi-polynomial for the group of reflections #w.r.t. x[i]-x[i+1]=-1 and x[k]=-1, according to the #product formula . For example, try: #Phi3F(t,2); Phi3F:=proc(t,k) local i,j: mul(mul((t[j]-t[i])*(1-t[i]*t[j]),j=i+1..k),i=1..k)/ mul(t[i]^(2*k-i+1),i=1..k)* mul(t[i]^2-1,i=1..k) : end: ###End Gessel-Zeilberger Formula#### ####Begin generalized Bessel Functions #JatN(a,t,N): The first N terms of Sum(t^(2*n+a)/n!/(n+a)!, n=0..N); #For example, try: JatN(1,t,10); JatN:=proc(a,t,N) local n,gu: gu:=0: for n from 0 while 2*n+a<=N do gu:=gu+t^(2*n+a)/n!/(n+a)!: od: gu: end: #MonoToJ(M,t,k,J): Given a monomial M in t[1], ..., t[k] #t[1]^a1*...*t[k]^ak, and a symbol J outputs J[a1]*...J[ak] #For example, try: #MonoToJ(t[1]^2*t[2],t,2,J); MonoToJ:=proc(M,t,k,J) local lu,i,d,gu: gu:=M: lu:=1: for i from 1 to k do d:=degree(M,t[i]): lu:=lu*J[abs(d)]: gu:=normal(gu/t[i]^d): od: if not type(gu,integer) then RETURN(FAIL): fi: lu*gu: end: #PolyToJ(P,t,k,J): Given a polynomial P in t[1], ..., t[k] #and a symbol J outputs the umbral image t^a->J[a] #For example, try: #PolyToJ(t[1]+t[2],t,2,J); PolyToJ:=proc(P,t,k,J) local i,P1: P1:=expand(P): if not type(P1,`+`) then RETURN(MonoToJ(P1,t,k,J)): else add(MonoToJ(op(i,P1),t,k,J),i=1..nops(P1)) : fi: end: #### ####GuessHolo2 without ezra with(SolveTools): #####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 #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 probabilites #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 remeber: 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 representaion #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 representaion #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 commuators 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 paramter 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,MaxC))): ope2:=numer(normal(FindrecF(i->F(i,21),n,N,MaxC))): ope3:=numer(normal(FindrecF(i->F(i,22),n,N,MaxC))): 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,MaxC))): 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) annihilationg #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,i,j,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: #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 d1nops(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: #######################End GuessHolo2 without ezra #SIPUR(K,L,L1,n,N): all the stories for walks #in (a) Z^k, (b) x1>=0, ..., xk>=0 #(c) x1>=x2>= ...>=xk (d) x1>=x2>= ...>=xk>=0 #for k=1..K #For example, try: #SIPUR(4,50,5000,n,N): SIPUR:=proc(K,L,L1,n,N) local k: for k from 1 to K do Sipur0(k,L,L1,n,N): od: for k from 1 to K do Sipur1(k,L,L1,n,N): od: for k from 1 to K do Sipur2(k,L,L1,n,N): od: for k from 1 to K do Sipur3(k,L,L1,n,N): od: end: #POLYA(K): all the Polya constants for walks #in (a) Z^k, (b) x1>=0, ..., xk>=0 #(c) x1>=x2>= ...>=xk (d) x1>=x2>= ...>=xk>=0 #for k=1..K #For example, try: #POLYA(5): POLYA:=proc(K) local k,i,x: for k from 2 to K do print(`For walks in the (hyper-cubic) lattice of dimension`, k): print(`The prob. of returning to the origin walking in Z^k is: `): print( Pr0([0$k],[0$k]) ): print(): print(`The prob. of returning to the origin walking in Z^k staying in `): print([seq(x[i]>=0, i=1..k)], ` equals `): print( Pr1([0$k],[0$k]) ): print(): print(`The prob. of returning to the origin walking in Z^k staying in `): print([seq(x[i]>=x[i+1], i=1..k-1)], ` equals `): print( Pr2([0$k],[0$k]) ): print(): print(`The prob. of returning to the origin walking in Z^k staying in `): print([seq(x[i]>=x[i+1], i=1..k-1), x[k]>=0 ], ` equals `): print( Pr3([0$k],[0$k]) ): print(): od: end: #POLYAc(K): all the Polya constants for walks #in (a) Z^k, (b) x1>=0, ..., xk>=0 #(c) x1>=x2>= ...>=xk (d) x1>=x2>= ...>=xk>=0 #for k=2..K, in compact form #For example, try: #POLYAc(5): POLYAc:=proc(K) local k: [seq( [ Pr0([0$k],[0$k]),Pr1([0$k],[0$k]),Pr2([0$k],[0$k]),Pr3([0$k],[0$k])] , k=2..K)]: end: