##################################################################### ## Capone.txt Save this file as Capone.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read Capone.txt # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## DoronZeil at gmail dot com # ###################################################################### print(`First Written: April 2021: tested for Maple 2018 `): print(`Version: April 2, 2021: `): print(): print(`This is Capone.txt, A Maple package`): print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's article: `): print(`Some Deep and Original Questions about the "critical exponents" of Generalized Ballot Sequences `): print(): print(`The most current version is available on WWW at:`): print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Capone.txt .`): print(`Please report all bugs to: DoronZeil at gmail dot com .`): print(`--------------------`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(``): print(`--------------------`): print(``): print(`For a list of the supporting functions type: ezra1(); For help with a specific function type ezra(FunctionName); `): print(): print(`--------------------`): print(``): print(`For a list of the FindRec.txt functions type: ezraFindRec(); For help with a specific function type ezraFindRec(FunctionName); `): print(`--------------------`): print(``): print(`--------------------`): print(``): print(`For a list of the STORY functions type: ezraS(); For help with a specific function type ezra(FunctionName); `): print(`--------------------`): print(``): ezra1:=proc() if args=NULL then print(`The SUPPORTING procedures are`): print(` DistN, MyAsy, TotW, W, SeqTotW1, SAW2, SAW2t, SAW2h Zinn `): else ezra(args): fi: end: ezraS:=proc() if args=NULL then print(`The Story procedures are`): print(` Paper2d, Paper2dT, Paper3d, Paper3dT, Paper4d `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` Capone.txt: A Maple package for studying Multi-dimensional Al Capone sequenes`): print(`The MAIN procedures are: Cap, MyAsyM , SeqTotW, Theo, TheoT `): print(``): elif nargs=1 and args[1]=Cap then print(`Cap(k,N): Let k be a list of positive integers of length d say. Let M=lcm(k[1],...,k[d]):`): print(`The first N terms, starting at n=1 in the number of ways of getting to [n*M/k[1],..., n*M/k[d]]`): print(`from the origin, always staying in k[1]*x[1]>=k[2]*x[2], k[2]*x[2]>=k[3]*x[3], ...k[d-1]*x[d-1]>=k[d]*x[d]. When k=[1$d] it is the`): print(`number of Young tableaux of the shape n^d. When d=2 it is R.C. Lyness' Al Capone sequence`): print(`It also returns the connective constant for the unrestricted walk`): print(`Try: `): print(`Cap([1,1,2],10);`): elif nargs=1 and args[1]=DistN then print(`DistN(k,N): inputs a vecor k of positive integers, and a positive integer N outputs the set`): print(`of points [a1,a2,...,ak] such that k[1]*a[1]>=k[2]*a[2]>=..>=k[d]*a[d] and the sum is N. Try`): print(`DistN([1,1],10);`): elif nargs=1 and args[1]=MyAsy then print(`MyAsy(L,n,k): tries to fit the sequence L to be of the form C*mu^n*n^theta*(1+c1/n+...+ck/n^k), where mu is unknown. Try:`): print(`MyAsy(Cap([1,1],100)[1],n,4 );`): elif nargs=1 and args[1]=MyAsyM then print(`MyAsyM(L,mu,n,k): tries to fit the sequence L to be of the form C*mu^n*n^theta*(1+c1/n+...+ck/n^k). Try:`): print(`MyAsyM(Cap([1,1],100),n,4 );`): elif nargs=1 and args[1]=Paper2d then print(`Paper2d(A,N,k1,n): inputs a positive integer A>=1 outputs an article`): print(`about the (estimated) asymptotics of all 2-dimensional Al Capone sequences with slope [m1,m2], m1>=m2 and relatively prime relative prime from 1 to A.`): print(`Try:`): print(`Paper2d(2,100,4,n);`): elif nargs=1 and args[1]=Paper2dT then print(`Paper2dT(A,N,k1,n): inputs a positive integer A>=1 outputs an article`): print(`about the (estimated) asymptotics of all sequenes enumerating walks with n steps in the region m[1]*x[1]-m[2]*x[2]>=0 `): print(`for all [m1,m2], m1>=m2 and relatively prime relative prime from 1 to A.`): print(`Try:`): print(`Paper2dT(2,100,4,n);`): elif nargs=1 and args[1]=Paper3d then print(`Paper3d(A,N,k1,n): inputs a positive integer A>=1 outputs an article`): print(`about the (estimated) asymptotics of all 3-dimensional Al Capone sequences with slope [m1,m2,m3], m1>=m2>=m3`): print(` and relatively prime relative prime from 1 to A.`): print(`Try:`): print(`Paper3d(2,100,4,n);`): elif nargs=1 and args[1]=Paper3dT then print(`Paper3dT(A,N,k1,n): inputs a positive integer A>=1 outputs an article`): print(`about the (estimated) asymptotics of all sequenes enumerating walks with n steps in the region m[1]*x[1]-m[2]*x[2]>=0, m[2]*x[2]-m[3]*x[3]>=0 `): print(`for all [m1,m2,m3], m1>=m2>=m3 and relatively prime relative prime from 1 to A.`): print(`Try:`): print(`Paper3dT(2,100,4,n);`): elif nargs=1 and args[1]=Paper4d then print(`Paper4d(A,N,k1,n): inputs a positive integer A>=1 outputs an article`): print(`about the (estimated) asymptotics of all 4-dimensional Al Capone sequences with slope [m1,m2,m3,m4], m1>=m2>=m3>=m4`): print(` and relatively prime relative prime from 1 to A.`): print(`Try:`): print(`Paper4d(2,40,4,n);`): elif nargs=1 and args[1]=SeqTotW then print(`SeqTotW(k,N): Inputs a list k of pos. integers k, of size d, say outputs the sequence of length N whose n-tn entry is thetotal number of lattice walks in`): print(`the d-dimensinal Manhattan lattice with N steps that stay in k[1]*x[1]>=k[2]*x[2]>=...>=k[d]*x[d]. Try`): print(`SeqTotW([1,1,1],20);`): elif nargs=1 and args[1]=SeqTotW1 then print(`SeqTotW1(k,N): Inputs a list k of pos. integers k, of size d, say outputs the list of length sum(k) whose`): print(`i-th item is the lhe list of length N`): print(`total number of lattice walks in of length sum(k)*n+i`): print(`the d-dimensinal Manhattan lattice with n steps that stay in k[1]*x[1]>=k[2]*x[2]>=...>=k[d]*x[d]. Try`): print(`SeqTotW1([1,1],10);`): elif nargs=1 and args[1]=SAW2 then print(`SAW2(): The first 79 terms of the sequence enumerating n-step self-avoiding walks in the 2D square lattice. Taken from the OEIS. Try:`): print(`SAW2();`): elif nargs=1 and args[1]=SAW2h then print(`SAW2h(): The first 106 terms of the sequence enumerating n-step self-avoiding walks in the 2D HEXAGONAL lattice. Taken from the OEIS. Try: `): print(`SAW2h();`): elif nargs=1 and args[1]=SAW2t then print(`SAW2t(): The first 40 terms of the sequence enumerating n-step self-avoiding walks in the 2D triangular lattice. Taken from the OEIS. Try: `): print(`SAW2t();`): elif nargs=1 and args[1]=Theo then print(`Theo(k,N,k1,n): Inputs a list k of positive integers of length d say, and a positive integer N>=20. Let M=lcm(k[1],...,k[d]):`): print(`It returns a Theorem about`): print(`the first N terms, starting at n=1 in the number of ways of getting to [n*M/k[1],..., n*M/k[d]]`): print(`from the origin, always staying in k[1]*x[1]>=k[2]*x[2], k[2]*x[2]>=k[3]*x[3], ...k[d-1]*x[d-1]>=k[d]*x[d]. When k=[1$d] it is the`): print(`number of Young tableaux of the shape n^d. When d=2 it is R.C. Lyness' Al Capone sequence`): print(`It also returns the connective constant for the unrestricted walk, and the estimated critical exponent`): print(`Try:`): print(`Theo([1,1,1],100,4,n);`): elif nargs=1 and args[1]=TheoT then print(`TheoT(k,N,k1,n): Let k be a list of positive integers of length d say. Let M=lcm(k[1],...,k[d]):`): print(`It returns a Theorem about`): print(`the first N terms, starting at n=1 in the number of walking n steps`): print(`from the origin, always staying in k[1]*x[1]>=k[2]*x[2], k[2]*x[2]>=k[3]*x[3], ...k[d-1]*x[d-1]>=k[d]*x[d]. When k=[1$d] it is the`): print(`number of Young tableaux with at most d rows with n cells. The case d=2 is binomial(n,n/2) and the case d=3 is`): print(`Amitai Regev's result that it is given by the Motzkin numbers.`): print(`It also returns the setimates of the connective constants for the subesequences with n mod sum(k).`): print(`Try:`): print(`TheoT([1,2],100,4,n);`): print(`TheoT([1,1,1],100,4,n);`): print(`TheoT([2,1,1],100,4,n);`): elif nargs=1 and args[1]=TotW then print(`TotW(k,N): Inputs a list k of pos. integers k, of size d, say outputs the total number of lattice walks in`): print(`the d-dimensinal Manhattan lattice with N steps that stay in k[1]*x[1]>=k[2]*x[2]>=...>=k[d]*x[d]. Try`): print(`TotW([1,1],10);`): elif nargs=1 and args[1]=W then print(`W(a,x,F): inputs a list a of non-negative integers a of length d, say, a letter x, and a list of expressions`): print(`in x[1],..., x[d], outputs the number of walks from the origin to the point a in the d-dim hyper-cubic lattice`): print(`where the expresions in F are all non-negative.. Try:`): print(`W([2,2,2],x,{x[1]-x[2],x[2]-x[3]});`): elif nargs=1 and args[1]=Zinn then print(`Zinn(L): inputs an integer sequence,uses Zinn-Justin's estimate for the connective constant and critical exponent.`): else print(`There is no such thing as`, args): fi: end: #DATA FROM THE OEIS about self-avoiding walks in 2D #NUMBER OF SELF-AVOIDING WALKS in the 2D square lattice SAW2:=proc(): [4, 12, 36, 100, 284, 780, 2172, 5916, 16268, 44100, 120292, 324932, 881500, 2374444, 6416596, 17245332, 46466676, 124658732, 335116620, 897697164, 2408806028, 6444560484,17266613812,46146397316,123481354908,329712786220,881317491628,2351378582244, 6279396229332,16741957935348,44673816630956,119034997913020,317406598267076,845279074648708,2252534077759844,5995740499124412, 15968852281708724,42486750758210044,113101676587853932,300798249248474268,800381032599158340,2127870238872271828,5659667057165209612, 15041631638016155884,39992704986620915140,106255762193816523332,282417882500511560972,750139547395987948108,1993185460468062845836, 5292794668724837206644,14059415980606050644844,37325046962536847970116,99121668912462180162908,263090298246050489804708, 698501700277581954674604,1853589151789474253830500,4920146075313000860596140,13053884641516572778155044,34642792634590824499672196, 91895836025056214634047716,243828023293849420839513468,646684752476890688940276172,1715538780705298093042635884,4549252727304405545665901684, 12066271136346725726547810652,31992427160420423715150496804,84841788997462209800131419244,224916973773967421352838735684,596373847126147985434982575724, 1580784678250571882017480243636,4190893020903935054619120005916,11107224538074654820152678182884,29442884996760677051402398150644, 78023796077779727644807609460228,206797849568186990141402577046860,547952781764285893561169365957068,1452142167241575828091155500636684, 3847327231644550282490410907667972,10194710293557466193787900071923676]: end: #NUMBER OF SELF-AVOIDING WALKS in the TRIANGULAR LATTICE SAW2t:=proc(): [1, 6,30,138,618,2730,11946,51882,224130,964134,4133166,17668938,75355206,320734686,1362791250,5781765582,24497330322,103673967882,438296739594, 1851231376374,7812439620678,32944292555934,138825972053046,584633909268402,2460608873366142,10350620543447034,43518414461742966, 182885110185537558,768238944740191374,3225816257263972170,13540031558144097474,56812878384768195282,238303459915216614558, 999260857527692075370,4188901721505679738374,17555021735786491637790,73551075748132902085986,308084020607224317094182,1290171266649477440877690, 5401678666643658402327390,22610911672575426510653226]: end: #NUMBER OF SELF-AVOIDING WALKS in the HEXAGONAL LATTICE SAW2h:=proc(): [1, 3, 6, 12, 24, 48, 90, 174, 336, 648, 1218, 2328, 4416, 8388, 15780, 29892, 56268, 106200, 199350, 375504, 704304, 1323996, 2479692, 4654464, 8710212, 16328220, 30526374, 57161568, 106794084, 199788408, 372996450, 697217994, 1300954248, 2430053136, 4531816950, 8459583678, 15769091448, 29419727280, 54816035922, 102216080286, 190380602052, 354843312276, 660671299170, 1230891734724, 2291023353264, 4266787588320, 7939282155480, 14780995214220, 27495750661500, 51174251554488, 95170708047492, 177077531456064, 329240750312604, 612431763064176, 1138447386722958, 2117139958203462, 3934744806805872, 7315657346364324, 13593701861615208, 25268621476344132, 46944923841500892, 87246030676633968, 162061613297503182, 301131372882048918, 559270924493111940, 1039015370238277392, 1929407176565844378, 3583868840338449552, 6654161548718435352, 12358172493655102908, 22942361044377497694, 42602546413521529998, 79079716539839010264, 146825880664057317624, 272508881502969583566, 505896263725614803028, 938837932546685757324, 1742680310186831144976, 3233705013725190208734, 6001731126079431880194, 11135634559782281744052, 20665346617856769531612, 38338839877458353531538, 71141110733855299445400, 131970447944873171324868, 244858119569459189454900, 454184324867092317646614, 842612982179925373445022, 1562820314973080810403960, 2899112977572410602989756, 5376639161689801476643794, 9973069503777017488350450, 18494428908117621975195672, 34302233761126805310473580, 63606602407525651570220172, 117963799130633513553560700, 218724567793338636147540660, 405611866822539347214154608, 752020561694426306189853768, 1394474897269109512017053416, 2585241775338665937474950376, 4793484596266980631401036304, 8886163328916481990779683742, 16475357899196133204173313054, 30540161675131744242797111332, 56619161516198581143528304882]: end: #END OF DATE FROM OEIS ##from FindRec ezraFindRec:=proc() if args=NULL then print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`): print(`equations satisfied by Discrete Functions of ONE Variable`): print(): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: `): print(` findrec, Findrec, FindrecF`): print(): elif nargs=1 and args[1]=findrec then print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`): print(`the sequence f of degree DEGREE and order ORDER.`): print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`): print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `): print(`of maximum DEGREE+ORDER<=MaxC`): print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`): elif nargs=1 and args[1]=FindrecF then print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`): print(`poly coffs. .g. try FindrecF(i->i!,n,N);`): elif nargs=1 and args[1]=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): print(`For example, try:`): print(`SeqFromRec(N-n-1,n,N,[1],10);`): fi: end: ###Findrec #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then RETURN(FAIL): fi: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 0 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #end Findrec with(linalg): #findrecEx(f,DEGREE,ORDER,m1): Explores whether thre #is a good chance that there is a recurrence of degree DEGREE #and order ORDER, using the prime m1 #For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003); findrecEx:=proc(f,DEGREE,ORDER,m1) local ope,var,eq,i,j,n0,eq1,a,A1, D1,E1,Eq,Var,f1,n,N: option remember: f1:=f mod m1: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1: od: eq:= eq union {eq1}: od: Eq:= convert(eq,list): Var:= convert(var,list): D1:=nops(Var): E1:=nops(Eq): if E10 then RETURN(false): fi: if E1-D1>=1 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+1],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: if E1-D1>=2 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+2],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: true: end: ##End from FindRec Digits:=100: #W(a,x,F): inputs a list a of non-negative integers a of length d, say, a letter x, and a list of expressions #in x[1],..., x[d], outputs the number of walks from the origin to the point a in the d-dim hyper-cubic lattice #where the expresions in F are all non-negative.. Try: #W([2,2,2],x,{x[1]-x[2],x[2]-x[3]}); W:=proc(a,x,F) local d,i,f: option remember: d:=nops(a): for i from 1 to d do if a[i]<0 then RETURN(0): fi: od: for f in F do if subs({seq(x[i]=a[i],i=1..d)},f)<0 then RETURN(0): fi: od: if a=[0$d] then RETURN(1): fi: add(W([op(1..i-1,a),a[i]-1, op(i+1..nops(a),a)],x,F),i=1..d): end: #Cap(k,N): Let k be a list of positive integers of length d say. Let M=lcm(k[1],...,k[d]): #The first N terms, starting at n=1 in the number of ways of getting to [n*M/k[1],..., n*M/k[d]] #from the origin, always staying in k[1]*x[1]>=k[2]*x[2], k[2]*x[2]>=k[3]*x[3], ...k[d-1]*x[d-1]>=k[d]*x[d]. When k=[1$d] it is the #number of Young tableaux of the shape n^d. When d=2 it is R.C. Lyness' Al Capone sequence #It also returns the connective constant for the unrestricted walk #Try: #Cap(3,[1,1,2],10); Cap:=proc(k,N) local d,i,F,gu,a,x,n,M,k1,mu: d:=nops(k): F:={seq(k[i]*x[i]-k[i+1]*x[i+1],i=1..d-1)}: M:=lcm(op(k)): k1:=[seq(M/k[i],i=1..nops(k))]: gu:=[]: for n from 1 to N do a:=[seq(k1[i]*n,i=1..d)]: gu:=[op(gu),W(a,x,F)]: od: mu:=convert(k1,`+`)^convert(k1,`+`)/mul(k1[i]^k1[i],i=1..d): gu,mu: end: #MyAsy(L,n,k): tries to fit the sequence L to be of the form C*mu^n*n^theta*(1+c1/n+...+ck/n^k). Try: MyAsy:=proc(L,n,k) local gu,logC,logmu,theta,c,i,eq,var,C,mu,x,ku: if nops(L)<20 then print(`List must be at least of length 20`): RETURN(FAIL): fi: if 10+k>nops(L) then print(`k must be at most`, nops(L)-10 ): fi: gu:=logC+logmu*n+theta*log(n)+add(c[i]/n^i,i=1..k): var:={logC,logmu,theta,seq(c[i],i=1..k)}: eq:={seq(log(L[i])-subs(n=i,gu),i=nops(L)-k-2..nops(L))}: var:=evalf(solve(eq,var)): C:=exp(subs(var,logC)): mu:=exp(subs(var,logmu)): theta:=subs(var,theta): C:=evalf(C,10): mu:=evalf(mu,10): theta:=evalf(theta,10): ku:= exp(subs(var,add(c[i]*x^i,i=1..k))): ku:=taylor(ku,x=0,k+1): ku:=add(evalf(coeff(ku,x,i),10)/n^i,i=0..k): evalf(C*mu^n*n^theta*ku,10): end: #MyAsyM(L,mu,n,k): tries to fit the sequence L to be of the form C*mu^n*n^theta*(1+c1/n+...+ck/n^k), where mu is given. Try: #gu:=Cap([2,3],200): MyAsyM(gu,n,3); MyAsyM:=proc(L,mu,n,k) local gu,logC,logmu,theta,c,i,eq,var,C,x,ku: if nops(L)<20 then print(`List must be at least of length 20`): RETURN(FAIL): fi: if 10+k>nops(L) then print(`k must be at most`, nops(L)-10 ): fi: logmu:=evalf(log(mu)): gu:=logC+logmu*n+theta*log(n)+add(c[i]/n^i,i=1..k): var:={logC,theta,seq(c[i],i=1..k)}: eq:={seq(log(L[i])-subs(n=i,gu),i=nops(L)-k-1..nops(L))}: var:=evalf(solve(eq,var)): C:=exp(subs(var,logC)): theta:=subs(var,theta): C:=evalf(C,10): theta:=evalf(theta,10): ku:= exp(subs(var,add(c[i]*x^i,i=1..k))): ku:=taylor(ku,x=0,k+1): ku:=add(evalf(coeff(ku,x,i),10)/n^i,i=0..k): evalf(theta,10),evalf(C,10)*mu^n*n^evalf(theta,10)*evalf(ku,10): end: #Theo(k,N,k1,n): Let k be a list of positive integers of length d say. Let M=lcm(k[1],...,k[d]): #It returns a Theorem about #the first N terms, starting at n=1 in the number of ways of getting to [n*M/k[1],..., n*M/k[d]] #from the origin, always staying in k[1]*x[1]>=k[2]*x[2], k[2]*x[2]>=k[3]*x[3], ...k[d-1]*x[d-1]>=k[d]*x[d]. When k=[1$d] it is the #number of Young tableaux of the shape n^d. When d=2 it is R.C. Lyness' Al Capone sequence #It also returns the connective constant for the unrestricted walk, and the estimated critical exponent #Try: #Theo([1,1,2],100,4,n); Theo:=proc(k,N,k1,n) local kc,d,gu,mu,theta,M,i,x,asy1: d:=nops(k): if not ( type(k,list) and type(N,integer) and type(N,integer) and N>=20 and min(op(k))>=1) then RETURN(FAIL): fi: M:=lcm(op(k)): kc:=[seq(M/k[i],i=1..nops(k))]: gu:=Cap(k,N): mu:=gu[2]: asy1:=MyAsyM(gu,n,k1): theta:=asy1[1]: asy1:=asy1[2]: gu:=gu[1]: print(`Theorem : Let a(n) be the number of ways of walking, in the`, nops(k) , `-dimensional Manhattan lattice with unit positive steps from `, [0$d], `to the point `): print(``): print([seq(kc[i]*n,i=1..d)]): print(``): print(`Always staying in the region `): print(``): print(seq( k[i]*x[i]-k[i+1]*x[i+1]>=0,i=1..d-1)): print(``): print(`The first`, N , `terms of this sequence are`): print(``): print(gu): print(``): print(`The connective constant is`, mu ): print(``): print(`The estimated critical expoent is`, theta ): print(``): print(`The estimated asymptotics is`): print(``): lprint(asy1): print(``): theta: end: #Paper2d(A,N,k1,n): inputs a positive integer A>=1 outputs an article #about the (estimated) asymptotics of all 2-dimensional Al Capone sequences with slope [k1,k2], k1>=k2 and relatively prime relative prime from 1 to A#Try: #Paper2d(2,100,4,n); Paper2d:=proc(A,N,k1,n) local gu,m1,m2,t0: t0:=time(): gu:=[]: for m1 from 1 to A do for m2 from 1 to m1 do if igcd(m1,m2)=1 then gu:=[op(gu), [[m1,m2],Theo([m1,m2],N,k1,n)]]: fi: od: od: print(`Here is a table of the critical exponents`): print(``): lprint(gu): print(``): print(`-------------------------------`): print(``): print(` This took `, time()-t0, `seconds. `): gu: end: #Paper2dT(A,N,k1,n): inputs a positive integer A>=1 outputs an article #about the (estimated) asymptotics of all 2-dimensional walks with n steps in the region k[1]*x-k[2]*y>=0 #Paper2dT(2,100,4,n); Paper2dT:=proc(A,N,k1,n) local gu,m1,m2,t0: t0:=time(): gu:=[]: for m1 from 1 to A do for m2 from 1 to m1 do if igcd(m1,m2)=1 then gu:=[op(gu), [[m1,m2],TheoT([m1,m2],N,k1,n)]]: fi: od: od: print(`Here is a table of the critical exponents`): print(``): lprint(gu): print(``): print(`-------------------------------`): print(``): print(` This took `, time()-t0, `seconds. `): gu: end: #Paper3d(A,N,k1,n): inputs a positive integer A>=1 outputs an article #about the (estimated) asymptotics of all 3-dimensional Al Capone sequences with slope [m1,m2,m3], m1>=m2>=m3 and relatively prime relative prime from 1 to A #Paper3d(2,100,4,n); Paper3d:=proc(A,N,k1,n) local gu,m1,m2,m3,t0: t0:=time(): gu:=[]: for m1 from 1 to A do for m2 from 1 to m1 do for m3 from 1 to m2 do if igcd(m1,m2,m3)=1 then gu:=[op(gu), [[m1,m2,m3],Theo([m1,m2,m3],N,k1,n)]]: fi: od: od: od: print(`Here is a table of the critical exponents`): print(``): lprint(gu): print(``): print(`-------------------------------`): print(``): print(` This took `, time()-t0, `seconds. `): gu: end: #Paper4d(A,N,k1,n): inputs a positive integer A>=1 outputs an article #about the (estimated) asymptotics of all 4-dimensional Al Capone sequences with slope [m1,m2,m3,m4], m1>=m2>=m3>=m4 #and relatively prime relative prime from 1 to A #Paper4d(2,100,4,n); Paper4d:=proc(A,N,k1,n) local gu,m1,m2,m3,m4,t0: t0:=time(): gu:=[]: for m1 from 1 to A do for m2 from 1 to m1 do for m3 from 1 to m2 do for m4 from 1 to m2 do if igcd(m1,m2,m3,m4)=1 then gu:=[op(gu), [[m1,m2,m3,m4],Theo([m1,m2,m3,m4],N,k1,n)]]: fi: od: od: od: od: print(`Here is a table of the critical exponents`): print(``): lprint(gu): print(``): print(`-------------------------------`): print(``): print(` This took `, time()-t0, `seconds. `): gu: end: #DistN(k,N): inputs a vecor k of positive integers, and a positive integer N outputs the set #of points [a1,a2,...,ak] such that k[1]*a[1]>=k[2]*a[2]>=..>=k[d]*a[d] and the sum is N DistN:=proc(k,N) local d,gu,k1,a1,lu,lu1: d:=nops(k): if N<0 then RETURN({}): fi: if N=0 then RETURN({[0$d]}): fi: if d=1 then RETURN({[N]}): fi: gu:={}: k1:=[op(2..d,k)]: for a1 from 0 to N do lu:=DistN(k1,N-a1): for lu1 in lu do if k[1]*a1>=k[2]*lu1[1] then gu:=gu union {[a1,op(lu1)]}: fi: od: od: gu: end: #TotW(k,N): Inputs a list k of pos. integers k, of size d, say outputs the total number of lattice walks in #the d-dimensinal Manhattan lattice with N steps that stay in k[1]*x[1]>=k[2]*x[2]>=...>=k[d]*x[d]. Try #TotW([1,1],10); TotW:=proc(k,N) local d,gu,a,x,F,i: gu:=DistN(k,N): d:=nops(k): F:={seq(k[i]*x[i]-k[i+1]*x[i+1],i=1..d-1)}: add(W(a,x,F), a in gu): end: #SeqTotW(k,N): Inputs a list k of pos. integers k, of size d, say outputs the list of length N whose n-th entry is #total number of lattice walks in #the d-dimensinal Manhattan lattice with n steps that stay in k[1]*x[1]>=k[2]*x[2]>=...>=k[d]*x[d]. Try #SeqTotW([1,1],10); SeqTotW:=proc(k,N) local n: [seq(TotW(k,n) ,n=1..N)]: end: #SeqTotW1(k,N): Inputs a list k of pos. integers k, of size d, say outputs the list of length sum(k) whose #i-th item is the lhe list of length N #total number of lattice walks in of length sum(k)*n+i #the d-dimensinal Manhattan lattice with n steps that stay in k[1]*x[1]>=k[2]*x[2]>=...>=k[d]*x[d]. Try #SeqTotW1([1,1],10); SeqTotW1:=proc(k,N) local gu,su,n,i: su:=convert(k,`+`): gu:=SeqTotW(k,su*(N+1)): [seq([seq(gu[su*n+i],n=0..N)],i=1..su)]: end: sn:=proc(resh,n1): -1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2): end: #Zinn(resh): Zinn-Justin's method to estimate #the C,mu, and theta such that #resh[i] is appx. Const*mu^i*i^theta #For example, try: #Zinn([seq(5*i*2^i,i=1..30)]); Zinn:=proc(resh) local s1,s2,theta,mu,n1,i: if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then RETURN(FAIL): fi: n1:=nops(resh)-1: s1:=sn(resh,n1): s2:=sn(resh,n1-1): theta:=evalf(2*(s1+s2)/(s1-s2)^2): mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): evalf([theta,mu],10): end: #TheoT(k,N,k1,n): Let k be a list of positive integers of length d say. Let M=lcm(k[1],...,k[d]): #It returns a Theorem about #the first N terms, starting at n=1 in the number of walking n steps #from the origin, always staying in k[1]*x[1]>=k[2]*x[2], k[2]*x[2]>=k[3]*x[3], ...k[d-1]*x[d-1]>=k[d]*x[d]. When k=[1$d] it is the #number of Young tableaux with at most d rows with n cells. The case d=2 is binomial(n,n/2) and the case d=3 is #Amitai Regev's result that it is given by the Motzkin numbers. #It also returns the setimates of the connective constants for the subesequences with n mod sum(k). #Try: #TheoT([1,2],100,4,n); #TheoT([1,1,1],100,4,n); #TheoT([2,1,1],100,4,n); TheoT:=proc(k,N,k1,n) local d,M,gu,i,gu1,lu,luE,asy1,theta,j,x: d:=nops(k): if not ( type(k,list) and type(N,integer) and type(N,integer) and N>=20 and min(op(k))>=1) then RETURN(FAIL): fi: M:=convert(k,`+`): gu:=SeqTotW(k,N): lu:=[]: luE:=[]: for i from 1 to M do gu1:=[seq(gu[i+M*j],j=0..trunc((nops(gu)-i)/M))]: asy1:=MyAsyM(gu1,d^M,n,k1): theta:=asy1[1]: lu:=[op(lu),theta]: luE:=[op(luE),asy1[2]]: od: print(``): print(`Theorem : Let a(n) be the number of ways of walking n steps, in the`, nops(k) , `-dimensional Manhattan lattice with unit positive steps `): print(``): print(`Always staying in the region `): print(``): print(seq( k[i]*x[i]-k[i+1]*x[i+1]>=0,i=1..d-1)): print(``): print(`The first`, N , `terms of this sequence, starting at n=1 are`): print(``): print(gu): print(``): print(`Assuming that The connective constants for the number of steps with i mod`, M, `for i from 1 to M are is`, d^M , `as it should ` ): print(``): print(`The estimated critical expoents for each of these moduli are`, lu ): print(``): print(`The estimated asymptotics are`): print(``): lprint(luE): print(``): lu: end: #Paper2dT(A,N,k1,n): inputs a positive integer A>=1 outputs an article #about the (estimated) asymptotics of all 2-dimensional walks with n steps in the region k[1]*x-k[2]*y>=0 #Paper2dT(2,100,4,n); Paper2dT:=proc(A,N,k1,n) local gu,m1,m2,t0: t0:=time(): gu:=[]: for m1 from 1 to A do for m2 from 1 to m1 do if igcd(m1,m2)=1 then gu:=[op(gu), [[m1,m2],TheoT([m1,m2],N,k1,n)]]: fi: od: od: print(`Here is a table of the critical exponents`): print(``): lprint(gu): print(``): print(`-------------------------------`): print(``): print(` This took `, time()-t0, `seconds. `): gu: end: #Paper3dT(A,N,k1,n): inputs a positive integer A>=1 outputs an article #with ThoremT([m1,m2,m3]) # for m1>=m2>=m3 and relatively prime relative prime from 1 to A #Paper3dT(2,100,4,n); Paper3dT:=proc(A,N,k1,n) local gu,m1,m2,m3,t0: t0:=time(): gu:=[]: for m1 from 1 to A do for m2 from 1 to m1 do for m3 from 1 to m2 do if igcd(m1,m2,m3)=1 then gu:=[op(gu), [[m1,m2,m3],TheoT([m1,m2,m3],N,k1,n)]]: fi: od: od: od: print(`Here is a table of the critical exponents`): print(``): lprint(gu): print(``): print(`-------------------------------`): print(``): print(` This took `, time()-t0, `seconds. `): gu: end: