(* ::Package:: *)

(* ::Text:: *)
(*STADJE.m is a package written by Luis A. Medina and Doron Zeilberger.*)
(*It is one of the packages that accompany the article :*)
(*'An Experimental Mathematics Perspective on the Old, and still Open, Question of When To Stop?'*)
(*	by Luis A. Medina and Doron Zeilberger.*)
(**)
(* Please report bug to : lmedina at math dot rutgers dot edu.*)
(**)
(*Please save STADJE.m in a folder that is readable by Mathematica.  For example, *)
(*if you are a Windows user, you can save the file STADJE.m in the following location:  C:\ProgramData\Mathematica\Applications.*)
(* *)
(*The most current version of this  package and article*)
(*are  available from http://www.math.rutgers.edu/~lmedina.*)
(**)
(* *)
(* Created: June 20, 2009.*)


BeginPackage["STADJE`"];


Mensaje::usage=" Created: June 20, 2009.

STADJE.m was written by Luis A. Medina and Doron Zeilberger,
Rutgers University. [medina, zeilberg] at math dot rutgers dot edu 

This is one of the packages that accompany the article:
'An Experimental Mathematics Perspective on the Old, and 
still Open, Question of When To Stop?'
	by Luis A. Medina and Doron Zeilberger.

Please report bug to: lmedina at math dot rutgers dot edu.

The most current version of this  package and paper
are  available from http://www.math.rutgers.edu/~lmedina.

For a list of the procedures simply input ?Help to see them."


Walks::usage="Walks[x,y,a,b] is the number of walks with unit steps [1,0],[0,1]
staying in the region y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x.

For example, try: Walks[5,10,1,1].";


ProbEsc::usage="ProbEsc[a,b,K] is the probability of escaping the region y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x 
in less than (a+b)*K steps.

For example, try:  ProbEsc[5,3,100].";


ExpectedGain::usage="ExpectedGain[a,b,K] is 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: `ExpectedGain[5,3,100]."


NetGain::usage="NetGain[a,b,K] is the relative net gain of keep going if you a 
heads and b tails right now in less than (a+b)*K steps."


WalksS::usage="WalksS[x,y,a,b] is the number of walks with unit steps [1,0],[0,1] 
staying in the region y > \!\(\*FractionBox[\"b\", \"a\"]\)x.

For example, try: WalksS[5,10,1,1]."


ProbTouch::usage="ProbTouch[a,b,K] is the probability of achieveing region 
y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x in less than (a+b)*K steps.

For example, try ProbTouch[5,3,100]."


NProbTouch::usage="NProbTouch[a,b,\[Epsilon]] is the probability of passing or
touching the line y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x in floating point with the 
current accuracy.

For example, try NProbTouch[5,3,\!\(\*SuperscriptBox[\"10\", 
RowBox[{\"-\", \"11\"}]]\)]."


NProbEsc::usage="NProbEsc[a,b,\[Epsilon]] is the probability of escaping the 
region y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x in floating point with the 
current accuracy.

For example, try: NProbEsc[5,3,\!\(\*SuperscriptBox[\"10\", 
RowBox[{\"-\", \"11\"}]]\)]."


AllProb::usage="AllProb[K,\[Epsilon]] returns the sequences for Table[Walks[i,i,a,b],{i,1,K}] 
(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 K terms.

For example, try: AllProb[4,\!\(\*SuperscriptBox[\"10\", 
RowBox[{\"-\", \"11\"}]]\)]."


MiddleAB::usage="MiddleAB[a,b,K,n] returns the first K terms of the sequence Walks[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]."


LineAB::usage="LineAB[a,b,K,n] returns 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]."


AllMiddleAB::usage="AllMiddleAB[M,K,L,n] returns the sequences for Table[Walks[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.

For example, try: AllMiddleAB[3,100,20,n]."


AllLineAB::usage="AllLineAB[M,K,L,n] returns the sequences for Table[Walks[a*i,b*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.

For example, try: AllLineAB[3,100,20,n]."



Help::usage="The functions included in this package are:
ExpectedGain
MiddleAB
LineAB
ProbEsc 
NProbEsc
ProbTouch
NProbTouch
Walks
WalksS
AllLineAB
AllMiddleAB
AllProb

If you want to know more information about these functions, 
simpliy input ?nameoffunction.  For example, try ?ProbEsc."


Begin["Private`"];


(* ::Text:: *)
(*The number of walks with unit steps [1,0],[0,1]*)
(*staying in the region y>=b/a*x. *)


Walks[x_Integer,y_Integer,a_Integer,b_Integer]/;x<0:=0;
Walks[0,y_Integer,a_Integer,b_Integer]:=1;
Walks[x_Integer,y_Integer,a_Integer,b_Integer]/;(y<b/a x):=0;
Walks[x_Integer,y_Integer,a_Integer,b_Integer]:=Walks[x,y,a,b]=Block[{$RecursionLimit=Infinity},Walks[x-1,y,a,b]+Walks[x,y-1,a,b]]


(* ::Text:: *)
(*All the pairs [i1,j1] such that [a*i+i1,b*i+j1]->[a*i+i1+1,b*i+j1] is the exiting edge*)
(*in the gambler's escape from y>=b/a*x *)


GettingOut[a_,b_]:=Table[{Floor[(a j)/b],j},{j,0,b-1}]


(* ::Text:: *)
(*The probability of escaping the region y>=b/a*x in less than (a+b)*K steps*)


ProbEsc[a_,b_,K_]:=Module[{go,s=0,count=0,g},
 go=GettingOut[a,b];
 While[(count=count+1)<=Length[go],
 g=go[[count]];
 s=s+Sum[Walks[a i+First[g],b i+Last[g],a,b]/2^((a+b)i+First[g]+Last[g]+1),{i,0,K}]
 ];
s
]


(* ::Text:: *)
(*The probability of escaping the region y>=b/a*x *)
(*in less than (a+b)*K steps in floating point*)


NProbEsc1[a_,b_,K_]:=Module[{go,g,count=0,s=0.0},
 go=GettingOut[a,b];
 While[(count=count+1)<=Length[go],
 g=go[[count]];
 s=s+Sum[Walks[a i+First[g],b i+Last[g],a,b]/2^((a+b)i+First[g]+Last[g]+1),{i,0,K}]
 ];
s
]


(* ::Text:: *)
(*The expected gain of keep going if you a heads and b tails right now in less than (a+b)*K steps*)


ExpectedGain[a1_Integer,b1_Integer,K_Integer]/;(a1>=0&&b1>=0&&K>=0):=Module[{go,g,rat,a,b,count=0,s=0},
 rat=a1/b1;
 a=Numerator[rat];
 b=Denominator[rat];
 go=GettingOut[a,b];
 While[(count=count+1)<=Length[go],
 g=go[[count]];
 s=s+Sum[Walks[a i+First[g],b i +Last[g],a,b]/2^((a+b)i+First[g]+Last[g]+1) ((a1+a i+First[g]+1)/(a1+b1+(a+b) i+First[g]+Last[g]+1)),{i,0,K}]
 ];
s+(1-ProbEsc[a,b,K])/2
]


(* ::Text:: *)
(*The relative net gain of keep going if you a heads and b tails right now*)
(*in less than (a+b)*K steps*)


NetGain[a_Integer,b_Integer,K_Integer]/;(a>=0&&b>=0&&K>=0):=(ExpectedGain[a,b,K]-a/(a+b))/(a/(a+b))


(* ::Text:: *)
(*The probability sequences (summands) of escaping the region y>=b/a*x in less than (a+b)*K steps*)


ProbSeqs[a_,b_,K_]:=Module[{go,g,list={},count=0},
 go=GettingOut[a,b];
 While[(count=count+1)<=Length[go],
 g=go[[count]];
 list=Append[list,Table[Walks[a i+g[[1]],b i+g[[2]],a,b]/2^((a+b) i+g[[1]]+g[[2]]+1),{i,0,K}]]
 ];
list
]


(* ::Text:: *)
(*The strip with parameters alpha and L of length K emanating from [a,b] (a>=b)*)


Strip[a_,b_,\[Alpha]_,L_,K_]:=Module[{gu,R,R1},
 gu=Table[Ceiling[b/a i],{i,0,a-1}];
 R=Join[ConstantArray[b,a-1],Flatten[Table[Table[b*j+gu[[i]],{i,1,Length[gu]}],{j,1,Floor[K/Length[gu]]}]]];
 R1=Table[Floor[\[Alpha] R[[i]]]+L,{i,1,Length[R]}];
 Thread[{R,R1}]
]


(* ::Text:: *)
(*The number of walks with unit steps [1,0],[0,1] starting at a0,b0 and staying in the region R.*)


gWalks[a0_,b0_,x_,y_,R_]/;(x<a0||y<b0):=0;
gWalks[a0_,b0_,a0_,b0_,R_]:=1;
gWalks[a0_,b0_,x_,y_,R_]/;(x>Length[R]||y<R[[x]][[1]]||y>R[[x]][[2]]):=0;
gWalks[a0_,b0_,x_,y_,R_]:=gWalks[a0,b0,x,y,R]=Block[{$RecursionLimit=Infinity},
 gWalks[a0,b0,x-1,y,R]+gWalks[a0,b0,x,y-1,R]
];


(* ::Text:: *)
(*The boundary of the region R expressed as a set of pairs [ptInRegion,ptNotInRegion]*)


BoundaryOfRegion[R_]:=Module[{boundary={},L,i=0,j=0},
 While[(i=i+1)<=Length[R]-1,
  If[First[R[[i]]]<First[R[[i+1]]],
   boundary=Union[boundary,{{{i,First[R[[i]]]},{i+1,First[R[[i]]]}}}]
  ];
 ];
 L=Length[R];
 boundary=Union[boundary,Table[{{L,j},{L+1,j}},{j,First[R[[L]]],R[[L]][[2]]}]];
 While[(j=j+1)<=Length[R],
  boundary=Union[boundary,{{{j,R[[j]][[2]]},{j,R[[j]][[2]]+1}}}]
 ];
 boundary
]


(* ::Text:: *)
(*The probability of escaping the region R starting at a0,b0*)


gProbEsc[a0_,b0_,R_]/;(b0<R[[a0]][[1]]||b0>R[[a0]][[2]]):=0;
gProbEsc[a0_,b0_,R_]:=Module[{B,b,gu=0,c=0},
 B=BoundaryOfRegion[R];
 While[(c=c+1)<=Length[B],
  b=B[[c]];
  gu=gu+gWalks[a0,b0,First[b[[1]]],Last[b[[1]]],R]/2^(b[[2]][[1]]+b[[2]][[2]]-a0-b0)
 ];
 gu
]


(* ::Text:: *)
(*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*)


gExpectedGain[a_,b_,R_]/;(b<R[[a]][[1]] || b>R[[a]][[2]]):=$Failed;
gExpectedGain[a0_,b0_,R_]:=Module[{B,b,gu=0,c=0},
 B=BoundaryOfRegion[R];
 While[(c=c+1)<=Length[B],
  b=B[[c]];
  gu=gu+gWalks[a0,b0,First[b[[1]]],Last[b[[1]]],R]/2^(b[[2]][[1]]+b[[2]][[2]]-a0-b0) Max[1/2,b[[2]][[1]]/(b[[2]][[1]]+b[[2]][[2]])]
 ];
gu
]


(* ::Text:: *)
(*The region y>=b/a*x*)


FindRegion[a_,b_,K_]:=Table[{Ceiling[b/a i],K},{i,1,K}]


(* ::Text:: *)
(*The number of walks with unit steps [1,0],[0,1] staying in the region y>b/a*x.*)


WalksS[x_Integer,y_Integer,a_Integer,b_Integer]/;x<0:=0;
WalksS[x_Integer,y_Integer,a_Integer,b_Integer]/;x==0:=1;
WalksS[x_Integer,y_Integer,a_Integer,b_Integer]/;y<=b/a x:=0;
WalksS[x_Integer,y_Integer,a_Integer,b_Integer]:=WalksS[x,y,a,b]=Block[{$RecursionLimit=Infinity},
 WalksS[x-1,y,a,b]+WalksS[x,y-1,a,b]
]


(* ::Text:: *)
(*The probability of achieveing region y>=b/a*x in less than (a+b)*K steps.*)


ProbTouch[a_Integer,b_Integer,K_Integer]/;(a>=0&&b>=0&&K>=0):=Module[{gu,g,c=0,s=0},
 gu=GettingOut[a,b];
 gu=Rest[gu];
 While[(c=c+1)<=Length[gu],
  g=gu[[c]];
  s=s+Sum[WalksS[a i+g[[1]],b i+g[[2]],a,b]/2^((a+b)i+g[[1]]+g[[2]]+1),{i,0,K}]
 ];
 s+Sum[WalksS[a i-1,b i,a,b]/2^((a+b)i),{i,1,K}]+1/2
]


(* ::Text:: *)
(*The probability of passing or touching the line y=b/a*x , in Floating point*)


NProbTouch1[a_,b_,K_]:=Module[{gu,g,c=0,s=0},
gu=GettingOut[a,b];
gu=Rest[gu];
While[(c=c+1)<=Length[gu],
 g=gu[[c]];
 s=s+Sum[WalksS[a i+g[[1]],b i+g[[2]],a,b]/2.0^((a+b)i+g[[1]]+g[[2]]+1),{i,0,K}]
];
s+Sum[WalksS[a i-1,b i,a,b]/2.0^((a+b)i),{i,1,K}]+1/2
]


(* ::Text:: *)
(*The probability of escaping the region y>=b/a*x  in floating point with the current accuracy*)


NProbEsc[a_Integer,b_Integer,\[Epsilon]_]/;(a>=0&&b>=0):=Module[{P1,P2,K=1},
P1=ProbEsc[a,b,100];
P2=ProbEsc[a,b,110];
While[(K=K+1)<=20 && Abs[P1-P2]>\[Epsilon],
P1=P2;
P2=ProbEsc[a,b,100+20K]
];
If[K==21,$Failed,N[P2,10]]
]


(* ::Text:: *)
(*The probability of passing or touching the line y=b/a*x in floating point with the current accuracy*)


NProbTouch[a_Integer,b_Integer,\[Epsilon]_]/;(a>=0&&b>=0):=Module[{K=1,P1,P2},
P1=ProbTouch[a,b,100];
P2=ProbTouch[a,b,110];
While[(K=K+1)<=20 &&Abs[P1-P2]>\[Epsilon],
P1=P2;
P2=ProbTouch[a,b,100+20 K]
];
If[K==21,$Failed,N[P2,10]]
]


(* ::Text:: *)
(*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.*)


AllProb[K_Integer,\[Epsilon]_]/;K>=0:=Module[{a=1,b,t0=0,te,tt,PE,PT},
Print["This gives all the probabilities of"];
Print["  (i) escaping the region y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x"];
Print["  (ii) escaping the region y > \!\(\*FractionBox[\"b\", \"a\"]\)x"];
Print["with unit steps to the right and up in the discrete lattice"];
Print["each with probability \!\(\*FractionBox[\"1\", \"2\"]\) starting at the origin for for 2 \[LessEqual] b \[LessEqual] a \[LessEqual] ",K];
While[(a=a+1)<=K,
 b=0;
 While[(b=b+1)<=a,
  If[GCD[a,b]==1,
  {te,PE}=Timing[NProbEsc[a,b,\[Epsilon]]];
  {tt,PT}=Timing[NProbTouch[a,b,\[Epsilon]]];
  t0=t0+te+tt;
  Print["a=", a, ", b=", b,", (i)=", PE,", (ii)=", PT]
  ]
 ]
];
Print["This took ", t0, " seconds."]
]


(* ::Text:: *)
(*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.*)


MiddleAB[a_Integer,b_Integer,K_Integer,n_Symbol]/;(a>=0&&b>=0&&K>=0):=Module[{seq1,seq2,err,D1,C},
seq1=Table[Walks[i,i,a,b],{i,1,K}];
seq2=Table[N[(seq1[[i]] Sqrt[i])/4^i],{i,1,Length[seq1]}];
err=Abs[seq2[[Floor[K/2]]]-seq2[[K]]];
D1=Floor[Log[1/err]/Log[10]]-1;
If[D1>2,
C=N[Floor[10^(D1+1) Last[seq2]]/10^(D1+1),D1];
{seq1,C 4^n/Sqrt[n]},
{seq1,$Failed}
]
]


(* ::Text:: *)
(*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.*)


LineAB[a_Integer,b_Integer,K_Integer,n_]/;(a>=0&&b>=0&&K>=0):=Module[{seq1,seq2,err,D1,C,c},
seq1=Table[Walks[a i,b i,a,b],{i,1,K}];
c=(a+b)^(a+b)/(a^a b^b);
seq2=Table[N[seq1[[i]]/c^i i^(3/2)],{i,1,Length[seq1]}];
err=Abs[seq2[[Floor[K/2]]]-seq2[[K]]];
D1=Floor[Log[1/err]/Log[10.]]-1;
If[D1>1,
C=N[Floor[10^(D1+1) Last[seq2]]/10^(D1+1),D1];
{seq1,C c^n/n^(3/2)},
{seq1,$Failed}
]
]


(* ::Text:: *)
(*The sequences for [Walks(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.*)


AllMiddleAB[M_Integer,K_Integer,L_Integer,n_Symbol]/;K>=L:=Module[{a=1,b,t,t0=0,list},
 Print["Below you can see the sequences enumerating the"];
 Print["number of ways of walking from the origin with"];
 Print["positive unit steps, in the square lattice to the"];
 Print["point (i,i), all while staying in the region y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x"];
 Print["followed by the conjectured asymptotics (if found)"];
 Print["for 2 \[LessEqual] b \[LessEqual] a \[LessEqual] ",K];
 While[(a=a+1)<=M,
  b=0;
  While[(b=b+1)<=a,
   If[GCD[a,b]==1,
   {t,list}=Timing[MiddleAB[a,b,K,n]];
   If[Last[list]=!=$Failed,
   Print["a=",a,", b=", b,", seq.:",First[list][[1;;L]]];
   Print["asymp.:", Last[list]],
   Print["a=",a,", b=", b,", seq.:",First[list][[1;;L]]];
   t0=t0+t;
   ]
  ]
 ]
];
Print["This took ",t0," seconds."]
]


(* ::Text:: *)
(*The sequences for Walks(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.*)


AllLineAB[M_Integer,K_Integer,L_Integer,n_Symbol]/;K>=L:=Module[{a=1,b,t,t0=0,list},
Print["Below you can see the sequences enumerating the"];
Print["number of ways of walking from the origin with"];
Print["positive unit steps, in the square lattice to the"];
Print["point (a i,b i), all while staying in the region y \[GreaterEqual] \!\(\*FractionBox[\"b\", \"a\"]\)x"];
Print["followed by the conjectured asymptotics (if found)"];
Print["for 2 \[LessEqual] b \[LessEqual] a \[LessEqual] ",K];
While[(a=a+1)<=M,
 b=0;
 While[(b=b+1)<=a,
  If[GCD[a,b]==1,
  {t,list}=Timing[LineAB[a,b,K,n]];
  t0=t0+t;
  If[Last[list]=!=$Failed,
   Print["a=",a,", b=", b,", seq.:",First[list][[1;;L]]];
   Print["asymp.:", Last[list]],
   Print["a=",a,", b=", b,", seq.:",First[list][[1;;L]]];
   ]
  ]
 ]
];
Print["This took ",t0," seconds."]
]


?Mensaje


End[];
EndPackage[];