###################################################################### ##STADJE: Save this file as STADJE # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read STADJE # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #First posted: June 23, 2009 print(`First posted: June 23, 2009`): print(` This is STADJE `): print(`It is one of the packages that accompany the article `): print(`"An Experimental Mathematics Perspective on the Old and still`): print(`open, question of when to stop?"`): print(`by Luis Medina and and Doron Zeilberger`): print(`available from Zeilberger's website`): 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/ .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are: BR,FindR, GettingOut `): print(`gExGain, gPrEsc, gW,`): print(` NetGain, PrSeqs, Strip, Zinn `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: `): print(` ExGain,MiddleAB, LineAB, `): print(` PrEsc, PrEscF, PrTouch, PrTouchF, `): print(` SipurLineAB, SipurMiddleAB, SipurPr , `): print(` W, Ws `): print(` `): elif nops([args])=1 and op(1,[args])=BR then print(`BR(R): the boundary of the region R expressed`): print(`as a set of pairs [ptInRegion,ptNotInRegion]`): print(`For example, try:`): print(`BR([[1,3]$3]);`): elif nops([args])=1 and op(1,[args])=ExGain then print(`ExGain(a,b,K): the expected gain of keep going`): print(`if you currently have`): print(`a heads and b tails right now`): print(`in less than (a+b)*K steps`): print(`and you stop as soon as you get a better score`): print(`or die at age (a+b)*K, where you are allowed to collect`): print(`1/2 in case your score is less than that.`): print(`For example, try:`): print(`ExGain(5,3,100);`): elif nops([args])=1 and op(1,[args])=FindR then print(`FindR(a,b,K): the region y>=b/a*x`): print(`For example, try:`): print(`FindR(5,3,10);`): elif nops([args])=1 and op(1,[args])=GettingOut then print(`GettingOut(a,b): all the pairs [i1,j1] such that [a*i+i1,b*i+j1]->[a*i+i1+1,b*i+j1] is the exiting edge`): print(`in the gambler's escape from y>=b/a*x `): print(`For example, try:`): print(`GettingOut(5,3);`): elif nops([args])=1 and op(1,[args])=gExGain then print(`gExGain(a0,b0,R): the expected gain of travelling in region R`): print(`starting at a0,b0 where the payoff is`): print(`max(1/2,x/(x+y)) as soon as you get out.`): print(`For example, try:`): print(`gExGain(1,1,[[1,5]$5]);`): elif nops([args])=1 and op(1,[args])=gPrEsc then print(`gPrEsc(a0,b0,R): the probability of escaping the region R`): print(`starting at a0,b0`): print(`For example, try:`): print(`gPrEsc(1,1,[[1,5]$5]);`): elif nops([args])=1 and op(1,[args])=gW then print(`gW(a0,b0,x,y,R): the number of walks with unit steps [1,0],[0,1]`): print(`starting at a0,b0 and`): print(`staying in the region R. For example, try:`): print(`gW(1,1,5,5,[[1,5]$5]);`): elif nops([args])=1 and op(1,[args])=LineAB then print(`LineAB(a,b,K,n): the first K terms of the sequence`): print(`W(a*i,b*i,a,b) (walks from the origin to (a*i,b*i)`): print(` staying in the region`): print(`y>=b/a*x), followed by the conjectured asympotics in n.`): print(`For example, try:`): print(`LineAB(3,2,100,n);`): elif nops([args])=1 and op(1,[args])=MiddleAB then print(`MiddleAB(a,b,K,n): the first K terms of the sequence`): print(`W(i,i,a,b) (walks from the origin to (i,i) staying in the region`): print(`y>=b/a*x), followed by the conjectured asympotics in n.`): print(`For example, try:`): print(`MiddleAB(3,2,100,n);`): elif nops([args])=1 and op(1,[args])=NetGain then print(`NetGain(a,b,K): the relative net gain of keep going`): print(`if you`): print(`a heads and b tails right now`): print(`in less than (a+b)*K steps`): print(`For example, try:`): print(`NetGain(5,3,100);`): elif nops([args])=1 and op(1,[args])=PrEsc then print(`PrEsc(a,b,K): the probability of escaping the region y>=b/a*x `): print(`in less than (a+b)*K steps.`): print(`For example, try:`): print(`PrEsc(5,3,100);`): elif nops([args])=1 and op(1,[args])=PrEscF then print(`PrEscF(a,b,eps): the probability of escaping the region y>=b/a*x `): print(`in floating point with error eps`): print(`For example, try:`): print(`PrEscF(5,3,10^(-11));`): elif nops([args])=1 and op(1,[args])=PrTouch then print(`PrTouch(a,b,K): the probability of escaping the region y>=b/a*x `): print(`or at least touching it, besides the origin`): print(`in less than (a+b)*K steps.`): print(`For example, try:`): print(`PrTouch(5,3,100);`): elif nops([args])=1 and op(1,[args])=PrTouchF then print(`PrTouchF(a,b,eps): the probability of escaping the region y>=b/a*x `): print(`or at least touching it, besides the origin`): print(`in floating point, with error =b/a*x `): print(`in less than (a+b)*K steps.`): print(`For example, try:`): print(`PrSeqs(5,3,100);`): elif nops([args])=1 and op(1,[args])=SipurMiddleAB then print(`SipurMiddleAB(M,K,L,n): the sequences for`): print(`[W(i,i,a,b),i=1..L] `): print(`(walks from the origin to (i,i) staying in the region`): print(`y>=b/a*x), followed by the conjectured asympotics in n`): print(`for starting at the origin for M>=a>=b>=1 and gcd(a,b)=1`): print(`it displays the first L terms. `): print(`SipurMiddleAB(3,100,20,n);`): elif nops([args])=1 and op(1,[args])=SipurPr then print(`SipurPr(K,eps): `): print(`all the probabilities of escaping the region y>=b/a*x `): print(`and or y>b/a*x (except at the beginning, of course)`): print(`with unit steps to the right and up in the discrete lattice`): print(`starting at the origin for K>=a>=b>=1 and gcd(a,b)=1`): print(`with error eps`): print(`For example, try:`): print(`SipurPr(4,10^(-11));`): elif nops([args])=1 and op(1,[args])=Strip then print(`Strip(a,b,alpha,L,K): the strip with parameters alpha and L`): print(`of length K`): print(`emanating from [a,b] (a>=b)`): print(`For example, try:`): print(`Strip(5,3,1,10,200);`): elif nops([args])=1 and op(1,[args])=W then print(`W(x,y,a,b): the number of walks with unit steps [1,0],[0,1]`): print(`staying in the region y>=b/a*x. For example, try:`): print(`W(5,10,1,1);`): elif nops([args])=1 and op(1,[args])=Ws then print(`Ws(x,y,a,b): the number of walks with unit steps [1,0],[0,1]`): print(`staying in the region y>b/a*x. For example, try:`): print(`Ws(5,10,1,1);`): elif nops([args])=1 and op(1,[args])=Zinn then print(` Zinn's method `): print(`Zinn(IntegerSequence): Given an increasing sequence a(n) of positive integers , expressed`): print(`in terms of a list, estimates the theta and mu such that `): print(` a(n) is asympt. to n^(theta)*mu^n. The output is `): print(` theta, mu `): else print(`There is no ezra for`,args): fi: end: #W(x,y,a,b): the number of walks with unit steps [1,0],[0,1] #staying in the region y>=b/a*x. For example, try: #W(5,10,1,1); W:=proc(x,y,a,b) option remember: if x<0 then 0: elif x=0 then 1: elif y[a*i+i1+1,b*i+j1] is the exiting edge #in the gambler's escape from y>=b/a*x #For example, try: #GettingOut(5,3); GettingOut:=proc(a,b) local j1: [seq([trunc(a*j1/b),j1],j1=0..b-1)]: end: #PrEsc(a,b,K): the probability of escaping the region y>=b/a*x #in less than (a+b)*K steps #For example, try: #PrEsc(5,3,100); PrEsc:=proc(a,b,K) local i,gu,g: gu:=GettingOut(a,b): add(add(W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K): end: #PrEscF1(a,b,K): the probability of escaping the region y>=b/a*x #in less than (a+b)*K steps in floating point #For example, try: #PrEscF1(5,3,100); PrEscF1:=proc(a,b,K) local i,gu,g: gu:=GettingOut(a,b): add(add(W(a*i+g[1],b*i+g[2],a,b)/2.^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K): end: #ExGain(a,b,K): the expected gain of keep going #if you #a heads and b tails right now #in less than (a+b)*K steps #For example, try: #ExGain(5,3,100); ExGain:=proc(a1,b1,K) local i,gu,g,shever,a,b: shever:=a1/b1: a:=numer(shever): b:=denom(shever): gu:=GettingOut(a,b): add(add( W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1)* (a1+a*i+g[1]+1)/(a1+b1+(a+b)*i+g[1]+g[2]+1), g in gu),i=0..K) +(1-PrEsc(a,b,K))/2: end: #ExGainOld(a,b,K): the expected gain of keep going #if you #a heads and b tails right now #in less than (a+b)*K steps #For example, try: #ExGainOld(5,3,100); ExGainOld:=proc(a1,b1,K) local i,gu,g,shever,a,b: shever:=a1/b1: a:=numer(shever): b:=denom(shever): gu:=GettingOut(a,b): add(add( W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1)* (a1+a*i+g[1]+1)/(a1+b1+(a+b)*i+g[1]+g[2]+1), g in gu),i=0..K): end: #NetGain(a,b,K): the relative net gain of keep going #if you a heads and b tails right now #in less than (a+b)*K steps #For example, try: (if it is neg. stop, if it is positive, stop) #NetGain(5,3,100); NetGain:=proc(a,b,K) local gu: gu:=ExGain(a,b,K): (gu-a/(a+b))/(a/(a+b)): end: #PrSeqs(a,b,K): the probability sequences (summands) #of escaping the region y>=b/a*x #in less than (a+b)*K steps #For example, try: #PrSeqs(5,3,100); PrSeqs:=proc(a,b,K) local i,gu,g: gu:=GettingOut(a,b): [seq( [seq(W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1),i=0..K)], g in gu)]: end: #Strip(a,b,alpha,L,K): the strip with parameters alpha and L #of length K #emanating from [a,b] (a>=b) #For example, try: #Strip(5,3,1,10,200); Strip:=proc(a,b,alpha,L,K) local gu,R,R1,i,j: gu:=[seq(ceil(b/a*i),i=0..a-1)]: R:=[b$(a-1),seq(seq(b*j+gu[i],i=1..nops(gu)),j=1..trunc(K/nops(gu)))]: R1:=[seq(trunc(alpha*R[i])+L,i=1..nops(R))]: [seq([R[i],R1[i]],i=1..nops(R))]: end: #gW(a0,b0,x,y,R): the number of walks with unit steps [1,0],[0,1] #starting at a0,b0 and #staying in the region R. For example, try: #gW(1,1,5,5,[[1,5]$5]); gW:=proc(a0,b0,x,y,R) option remember: if xnops(R) or yR[x][2] then RETURN(0): else RETURN(gW(a0,b0,x-1,y,R)+gW(a0,b0,x,y-1,R)): fi: end: #BR(R): the boundary of the region R expressed #as a set of pairs [ptInRegion,ptNotInRegion] #For example, try: #BR([[1,3]$3]); BR:=proc(R) local gu,k,j,i: gu:={}: for i from 1 to nops(R)-1 do if R[i][1]R[a0][2] then RETURN(FAIL): fi: B:=BR(R): gu:=0: for b in B do gu:=gu+gW(a0,b0,op(b[1]),R)/2^(b[2][1]+b[2][2]-a0-b0): od: gu: end: #gExGain(a0,b0,R): the expected gain of travelling in region R #starting at a0,b0 where the payoff is #max(1/2,x/(x+y)) as soon as you get out #For example, try: #gExGain(1,1,[[1,5]$5]); gExGain:=proc(a0,b0,R) local B,b,gu: if b0R[a0][2] then RETURN(FAIL): fi: B:=BR(R): gu:=0: for b in B do gu:=gu+gW(a0,b0,op(b[1]),R)/2^(b[2][1]+b[2][2]-a0-b0)* max(1/2,b[2][1]/(b[2][1]+b[2][2])): od: gu: end: #FindR(a,b,K): the region y>=b/a*x #For example, try: #FindR(5,3,10); FindR:=proc(a,b,K) local i: [seq([ceil(b*i/a),K],i=1..K)]: end: #PrTouch(a,b,K): the probability of achieveing region y>=b/a*x #in less than (a+b)*K steps #For example, try: #PrTouch(5,3,100); PrTouch:=proc(a,b,K) local i,gu,g: gu:=GettingOut(a,b): gu:=[op(2..nops(gu),gu)]: add(add(Ws(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K) + add(Ws(a*i-1,b*i,a,b)/2^((a+b)*i),i=1..K)+1/2: : end: #PrTouchF1(a,b,K): the probability of achieveing region y>=b/a*x #in less than (a+b)*K steps, in Floating point #For example, try: #PrTouchF1(5,3,100); PrTouchF1:=proc(a,b,K) local i,gu,g: gu:=GettingOut(a,b): gu:=[op(2..nops(gu),gu)]: add(add(Ws(a*i+g[1],b*i+g[2],a,b)/2.^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K) + add(Ws(a*i-1,b*i,a,b)/2.^((a+b)*i),i=1..K)+1/2: : end: #Ws(x,y,a,b): the number of walks with unit steps [1,0],[0,1] #staying in the region y>b/a*x. For example, try: #Ws(5,10,1,1); Ws:=proc(x,y,a,b) option remember: if x<0 then 0: elif x=0 then 1: elif y<=b/a*x then RETURN(0): else RETURN(Ws(x-1,y,a,b)+Ws(x,y-1,a,b)): fi: end: #PrEscF(a,b): the probability of escaping the region y>=b/a*x #in floating point with the current accuracy #For example, try: #PrEscF(5,3); PrEscF:=proc(a,b,eps) local K,gu1,gu2: gu1:=PrEsc(a,b,100): gu2:=PrEsc(a,b,110): for K from 2 to 20 while abs(gu1-gu2)>eps do gu1:=gu2: gu2:=PrEsc(a,b,100+20*K): od: if K=21 then FAIL: else evalf(gu2): fi: end: #PrTouchF(a,b): the probability of passing or touching the line y=b/a*x #in floating point with the current accuracy #For example, try: #PrTouchF(5,3); PrTouchF:=proc(a,b,eps) local K,gu1,gu2: gu1:=PrTouch(a,b,100): gu2:=PrTouch(a,b,110): for K from 2 to 20 while abs(gu1-gu2)>eps do gu1:=gu2: gu2:=PrTouch(a,b,100+20*K): od: if K=21 then FAIL: else evalf(gu2): fi: end: #SipurPr(K,eps): #all the probabilities of escaping the region y>=b/a*x #and or y>b/a*x (except at the beginning, of course) #with unit steps to the right and up in the discrete lattice #starting at the origin for K>=a>=b>=1 and gcd(a,b)=1 #with error eps #For example, try: #SipurPr(4,10^(-11)); SipurPr:=proc(K,eps) local a,b,t0: t0:=time(): print(`This is the story of all the probabilities of`): print(`(i) escaping the region y>=b/a*x `): print(`(ii) escaping y>b/a*x (except at the beginning, of course)`): print(`with unit steps to the right and up in the discrete lattice`): print(`each with prob. 1/2 `): print(`starting at the origin`): print(`for 2<=b<=a<=`, K): for a from 2 to K do for b from 1 to a do if gcd(a,b)=1 then print(`a=`,a,`b=`,b, `(i)=`, PrEscF(a,b,eps),`(ii)=`, PrTouchF(a,b,eps)): fi: od: od: print(`------------------------`): print(`This took`, time()-t0, `seconds. `): 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): end: #MiddleAB(a,b,K,n): the first K terms of the sequence #W(i,i,a,b) (walks from the origin to (i,i) staying in the region #y>=b/a*x), followed by the conjectured asympotics in n #For example, try: #MiddleAB(3,2,100,n); MiddleAB:=proc(a,b,K,n) local gu,i,gu1,err,D1,C: gu:=[seq(W(i,i,a,b),i=1..K)]: gu1:=[seq(evalf(gu[i]/4^i*sqrt(i)),i=1..nops(gu))]: err:=abs(gu1[trunc(K/2)]-gu1[K]): D1:=trunc(log(1/err)/log(10.))-1: if D1>2 then C:=evalf(gu1[nops(gu1)],D1): RETURN(gu,C*4^n/sqrt(n)): else RETURN(gu,FAIL): fi: end: #LineAB(a,b,K,n): the first K terms of the sequence #W(a*i,b*i,a,b) (walks from the origin to (i,i) staying in the region #y>=b/a*x), followed by the conjectured asympotics in n #For example, try: #LineAB(3,2,100,n); LineAB:=proc(a,b,K,n) local gu,i,gu1,err,D1,C,c: gu:=[seq(W(a*i,b*i,a,b),i=1..K)]: c:=(a+b)^(a+b)/a^a/b^b: gu1:=[seq(evalf(gu[i]/c^i*i^(3/2)),i=1..nops(gu))]: err:=abs(gu1[trunc(K/2)]-gu1[K]): D1:=trunc(log(1/err)/log(10.))-1: if D1>1 then C:=evalf(gu1[nops(gu1)],D1): RETURN(gu,C*c^n/n^(3/2)): else RETURN(gu,FAIL): fi: end: #SipurMiddleAB(M,K,L,n): the sequences for #[W(i,i,a,b),i=1..L] (walks from the origin to (i,i) staying in the region #y>=b/a*x), followed by the conjectured asympotics in n #for starting at the origin for M>=a>=b>=1 and gcd(a,b)=1 #it displays the first L terms #SipurMiddleAB(3,100,20,n); SipurMiddleAB:=proc(M,K,L,n) local a,b,t0,gu: t0:=time(): print(`Below are the sequences`): print(`enumerating the number of ways of walking from the origin`): print(`with positive unit steps, in the square lattice`): print(`to the point (i,i), all the while staying in the`): print(`the region y>=b/a*x `): print(`followed by the conjectured asymptotics (if found) `): print(`for 2<=b<=a<=`, M): for a from 2 to M do for b from 1 to a do if gcd(a,b)=1 then gu:=MiddleAB(a,b,K,n): if gu[2]<>FAIL then print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]): print(`asymp.:`, gu[2]): else print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]): fi: fi: od: od: print(`------------------------`): print(`This took`, time()-t0, `seconds. `): end: #SipurLineAB(M,K,L,n): the sequences for #[W(i,i,a,b),i=1..L] (walks from the origin to (i,i) staying in the region #y>=b/a*x), followed by the conjectured asympotics in n #for starting at the origin for M>=a>=b>=1 and gcd(a,b)=1 #it displays the first L terms #SipurLineAB(3,100,20,n); SipurLineAB:=proc(M,K,L,n) local a,b,t0,gu: t0:=time(): print(`Below are the sequences`): print(`enumerating the number of ways of walking from the origin`): print(`with positive unit steps, in the square lattice`): print(`to the point (a*i,b*i), all the while staying in the`): print(`the region y>=b/a*x `): print(`followed by the conjectured asymptotics (if found) `): print(`for 2<=b<=a<=`, M): for a from 2 to M do for b from 1 to a do if gcd(a,b)=1 then gu:=LineAB(a,b,K,n): if gu[2]<>FAIL then print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]): print(`asymp.:`, gu[2]): else print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]): fi: fi: od: od: print(`------------------------`): print(`This took`, time()-t0, `seconds. `): end: