#####################################################################
## MultiPropp.txt Save this file as   MultiPropp.txt to use it,      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `MultiPropp.txt` <Enter>                     #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################
with(plots):
print(`First Written: Aug. 2026: tested for Maple 2025 `):
print(`This Version: Aug. 6, 2026  `):
print():
print(`This is  MultiPropp.txt, a Maple package`):
print(`one of the two Maple packages accompanying  Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(`Estimating  Many Constants with a  Coin  `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/MultiPropp.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`For a list of the less important procedures`):
print(` type "ezra1();". For specific help type "ezra(procedure_name);" `):

print(`For a list of the Story functions type: ezraSt();`):

print(`For a list of the Simulations functions type: ezraSi();`):

ezraSi:=proc()
if args=NULL then
print(`The Simulations procedures are:  Roll, ROU,  Simu, WalkOut  `):
else
ezra(args):
fi:
end:

ezraSt:=proc()
if args=NULL then
print(`The Story procedures are:  Paper1 `):
else
ezra(args):

fi:
end:

ezra1:=proc()
if args=NULL then
print(`The SUPPORTING procedures are: Pars, ParsE. YF`):
else
ezra(args):
fi:
end:
ezra:=proc()
if args=NULL then

print(` MultiPropp.txt: A Maple package for  estimating constants by tossing a coin and using Symbolic computation to study them by travelling in many dimensions, not just two.`):
print(`The MAIN procedures are: AppExp0  `):
print(``):


elif nargs=1 and args[1]=AppExp0 then
print(`AppExp0(P,f,x,K):`):
print(`Inputs a probabiity distribution,P, (let k:=nops(P)) on {1,..,k} where P[i] is the prob. that the walker walks in the step [0$(i-1),1,0$(k-i)], `):
print(`an expression f in x[1], ..., x[k], does the walk until it leavs  {x[i]-x[i+1]>=0,i=1..k-1} `):
print(`Output:  the approximate expected value of f(pt) where pt is the point where it is outside for the first time limited to at most K steps`):
print(`followed by the expected number of steps. Try:`):
print(`AppExp0([1/3,1/3,1/3],x[3]/(x[1]+x[2]+x[3]),x,12):`):

elif nargs=1 and args[1]=Paper1 then
print(`Paper1(K,K1,K2): a paper about the expected value, both the exact value (via the truncated sum after K2 terms ) and the simulation values for the expected value`):
print(`of (x[i])/(x[1]+x[2]+...+x[k]) until getting for the first time outside x[1]>=x[2]>=...x[k]>=0 for i from 1 to k and k from 2 to K.`):
print(`It also compares it to the simulation estimate by doint it K1 times. Try:`):
print(`Paper1(5,1000,100);`):

elif nargs=1 and args[1]=Pars then
print(`Pars(n,k): The set of [a1,a2,..., ak] with a1>=a2>=..ak>=0 and a1+a2+...ak=n. Try:`):
print(`Pars(10,3);`):

elif nargs=1 and args[1]=ParsE then
print(`ParsE(n,k,i): The set of [a1,a2,..., ak] with a1>=a2>=..ak>=0 and a1+a2+...ak=n and a[i]=a[i+1]. Try:`):
print(`ParsE(10,3,1);`):


elif nargs=1 and args[1]=Roll then
print(`Roll(P): rolls a die with prob. dist. P. Try:`):
print(`add(x[Roll([1/6,1/3,1/2])],i=1..600);`):

elif nargs=1 and args[1]=ROU then
print(`ROU(P): Given a loaded die with prob. dist. P finds the roullete. Try:`):
print(`ROU([1/3,1/6,1/2]);`):


elif nargs=1 and args[1]=Simu then
print(`Simu(P,d,f,x,K,LI): Inputs a probabiity distribution,P, (let k:=nops(P)) on {1,..,k} where P[i] is the prob. that the walker walks in the step [0$(i-1),1,0$(k-i)], `):
print(`an expression f in x[1], ..., x[k], does the walk until it leavs  {x[i]-x[i+1]>=-d,i=1..k-1} for the first time, limited to LI steps .  `):
print(`Output:  the triple consisting of`):
print(`(i) the average value of f(destination) for those that made it by LI steps, hence giving an estimate for the expectated value of f at leaving`):
print(`(ii) the average number of steps until it went out`):
print(`(iii) the number of times it made it to the goal. Try:`):
print(`Simu([1/4$4],0,(x[1]+1)/(x[1]+x[2]+x[3]+x[4]),x,100,1000);`):

elif nargs=1 and args[1]=WalkOut then
print(`WalkOut(P,d,LI): Inputs a prob. P walks for at most LI steps, until for the first time it is out of x[i]-x[i+1]>=-d. If it doesn't make it in LI steps it returns FAIL. Try:`):
print(`WalkOut([1/3,1/3,1/3],0,1000);`):

elif nargs=1 and args[1]=YF then
print(`YF(a): The number of walks from the origin to a in (k:=nops(a)) x[1]>=x[2]>=...x[k]>=0 using the Young-Frobenius formula. Try:`):
print(`YF([5,3,2]);`):

else
 print(`There is no such thing as`, args):

fi:
end:




#ROU(P): Given a loaded die with prob. dist. P finds the roullete. Try:
#ROU([1/3,1/6,1/2]);
ROU:=proc(P) local k,i,M:
option remember:
k:=nops(P):
if not (convert(P,`+`)=1 and min(P)>=0 and max(P)<=1) then
 RETURN(FAIL):
 fi:
M:=lcm(seq(denom(P[i]),i=1..k)):
[seq(M*P[i],i=1..k)]:
end:

#Roll1(L): rolls a roullette L
Roll1:=proc(L) local i,n,r,j:
n:=convert(L,`+`):
r:=rand(1..n)():
for i from 1 while add(L[j],j=1..i)<r do od:
i:
end:

#Roll(P): rolls a die with prob. dist. P. Try:
#add(x[Roll([1/6,1/3,1/2]),i=1..600);
Roll:=proc(P) : Roll1(ROU(P)):end:





#WalkOut(P,d,LI): Inputs a prob. P walks for at most LI steps, until for the first time it is out of x[i]-x[i+1]>=-d. If it doesn't make it in LI steps it returns FAIL. Try:
#WalkOut([1/3,1/3,1/3],0,1000);
WalkOut:=proc(P,d,LI) local k,W,pt,j,i:
k:=nops(P):
pt:=[0$k]:
W:=[pt]:

while min(seq(pt[i]-pt[i+1],i=1..k-1))>=-d and nops(W)<=LI do
j:=Roll(P):
pt:=pt+[0$(j-1),1,0$(k-j)]:
W:=[op(W),pt]:
od:

if nops(W)=LI+1 then
FAIL:
else
W:
fi:

end:


#Simu(P,d,f,x,K,LI): Inputs a probabiity distribution,P, (let k:=nops(P)) on {1,..,k} where P[i] is the prob. that the walker walks in the step [0$(i-1),1,0$(k-i)], 
#an expression f in x[1], ..., x[k], does the walk until it leavs  {x[i]-x[i+1]>=-d,i=1..k-1} for the first time, limited to LI steps .  
#Output:  the triple consisting of
#(i) the average value of f(destination) for those that made it by LI steps, hence giving an estimate for the expectated value of f at leaving
#(ii) the average number of steps until it went out
#(iii) the number of times it made it to the goal.
#Try:
#Simu([1/4$4],0,(x[1]+1)/(x[1]+x[2]+x[3]+x[4]),x,100,1000);
Simu:=proc(P,d,f,x,K,LI) local i,W,co,len,val,i1,W1:

co:=0:
len:=0:
val:=0:

for i from 1 to K do
W:=WalkOut(P,d,LI):
W1:=W[nops(W)]:
if W<>FAIL then
co:=co+1:
len:=len+nops(W)-1:
val:=val+subs({seq(x[i1]=W1[i1],i1=1..nops(P))},f):
fi:
od:
if co=0 then
 RETURN(FAIL):
fi:
evalf([val/co,len/co,co]):
end:




#YF(a): The number of walks from the origin to a in (k:=nops(a)) x[1]>=x[2]>=...x[k]>=0 using the Young-Frobenius formula. Try:
#YF([5,3,2]);
YF:=proc(a) local k,i,j:
k:=nops(a):
mul(mul(a[i]-a[j]+j-i,j=i+1..k),i=1..k)*add(a[i],i=1..k)!/mul((a[i]+k-i)!,i=1..k):
end:

#Pars(n,k): The set of [a1,a2,..., ak] with a1>=a2>=..ak>=0 and a1+a2+...ak=n. Try:
#Pars(10,3);
Pars:=proc(n,k) local gu,lu,lu1:
option remember:
if n<0 then
 RETURN({}):
fi:
if n=0 then
 RETURN({[0$k]}):
fi:
if k=1 then
 RETURN({[n]}):
fi:
lu:=Pars(n,k-1):
gu:={seq([op(lu1),0], lu1 in lu)}:
lu:=Pars(n-k,k):
gu union {seq(lu1+[1$k],lu1 in lu)}:
end:

#ParsE(n,k,i): The partitions of n into k parts (allowing 0) where a[i]=a[i+1], of course 1<=i<k . Try:
#ParsE(8,3,1);
ParsE:=proc(n,k,i) local lu,lu1,gu:
option remember:
if not (i>=1 and i<k) then
 RETURN(FAIL):
fi:
gu:={}:

lu:=Pars(n,k):

for lu1 in lu do
 if lu1[i]=lu1[i+1] then
  gu:=gu union {lu1}:
 fi:
od:
gu:
end:



#AppExp0(P,f,x,K):
#Inputs a probabiity distribution,P, (let k:=nops(P)) on {1,..,k} where P[i] is the prob. that the walker walks in the step [0$(i-1),1,0$(k-i)], 
#an expression f in x[1], ..., x[k], does the walk until it leavs  {x[i]-x[i+1]>=0,i=1..k-1} 
#Output:  the approximate expected value of f(pt) where pt is the point where it is outside for the first time limited to at most K steps
#followed by the expected number of steps. Try:
#AppExp0([1/3,1/3,1/3],x[3]/(x[1]+x[2]+x[3]),x,12):
AppExp0:=proc(P,f,x,K) local k,gu,i,i1,lu,n,lu1:
k:=nops(P):

if not (type(P,list) and convert(P,`+`)=1 and max(P)<=1 and min(P)>=0) then
print(P, `should be a prob. distribution`):
RETURN(FAIL):
fi:

gu:=0:

for n from 0 to K do
for i from 1 to k-1 do
lu:=ParsE(n,k,i):

lu:={seq(lu1+[0$i,1,0$(k-i-1)], lu1 in lu)}:

for lu1 in lu do


gu:=gu+YF(lu1-[0$i,1,0$(k-i-1)])*mul(P[i1]^lu1[i1],i1=1..k)*subs({seq(x[i1]=lu1[i1],i1=1..k)},f):

od:
od:
od:

gu:
end:



#Paper1(K,K1,K2): a paper about the expected value, both the exact value (via the truncated sum at K2) and the simulation values for the expected value
#of (x[i])/(x[1]+x[2]+...+x[k]) until getting for the first time outside x[1]>=x[2]>=...x[k]>=0 for i from 1 to k and k from 2 to K.
#It also compares it to the simulation estimate by doint it K1 times. Try:
#Paper1(5,1000,100);
Paper1:=proc(K,K1,K2) local gu1,gu2,k,i,i1,P,f,x:
print(``):
print(`The Expected value of x[i]/(x[1]+....+x[k]) for i from 1 to k and  k from 2 to `, K, `via exact (truncated) summation and via simulation`, K1, `times `):
print(``):
print(`of the first time the fair walker goes out of x[1]>=x[2]>=...>=x[k]>=0, starting at the origin`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

for k from 2 to K do
P:=[(1/k)$k]:

print(`In `, k , ` dimensions `):
for i from 1 to k do

f:=x[i]/add(x[i1],i1=1..k):

print(`The expected values followed by the simulated approximations for`, f, ` are respecively`):
gu1:=evalf(AppExp0(P,f,x,K2)) :

gu2:=Simu(P,0,f,x,K1,1000)[1]:

print(gu1,gu2):
od:
od:
print(``):
print(`---------------------------------------------`):
print(``):
print(`This took`, time(), `seconds. `):
print(``):

end: