#C27.txt, April 30, 2020 Dr. Z's Experimental Mathematics Class Math 640 (RU)
Help:=proc(): print(`RW(d,K), PolCons(d,K,M) `): end:

#RW(d,K): the number of steps for a  simple random walk in d-dimensional space until it returns to the origin
#it outputs the number of steps or FAIL
RW:=proc(d,K) local pt,co,i,j,ra1,ra2:
pt:=[0$d]:
ra1:=rand(1..d):
ra2:=rand(0..1):

for co from 1 to K do
i:=ra1():
j:=ra2():

if j=0 then
 pt:=[op(1..i-1,pt),pt[i]-1,op(i+1..d,pt)]:
else
 pt:=[op(1..i-1,pt),pt[i]+1,op(i+1..d,pt)]:
fi:

if pt=[0$d] then
 RETURN(co):
fi:

od:

FAIL:

end:

#PolCons(d,K,M): estimating Polya's random walk constant in d dimension, pretending that K is infinity
#it also estimates the expected time of return
#doing it M times
#Try:
#PolCons(3,10000,1000);
PolCons:=proc(d,K,M) local co1,co2,i,W:
co1:=0:
co2:=0:
for i from 1 to M do
  W:=RW(d,K):
   if W<>FAIL then
     co1:=co1+1:
      co2:=co2+W:
   fi:
od:

evalf([co1/M,co2/M]):

end: