Help:=proc(): print(` LD(p,x) , GenM(A,B,p), Fatten1(V,S) , Fatten(v,S)`): 
print(`RO(M), MC(A,B,K,p) , Wt(M,p)`):
end:


#LD(p,x): Simulates one throw of a die whose prob. gen. function
#is p(x). For example LD((x+x^2+x^3+x^4+x^5+x^6)/6,x);
LD:=proc(p,x) local M,P,i,L,L2,ra,r:

P:=expand(p/subs(x=1,p)):

M:=lcm( seq( denom(coeff(P,x,i)), i=ldegree(P,x)..degree(P,x))):


L:=[]:
for i from ldegree(P,x) to degree(P,x) do
 if coeff(P,x,i)<>0 then
    L:=[op(L), [i,M*coeff(P,x,i)]]:
 fi:
od:

L2:=[seq(L[i][2],i=1..nops(L))]:
L2:=[seq( add(L2[j],j=1..i),i=1..nops(L) )]:

ra:=rand(1..M):
r:=ra():

for i from 1 to nops(L2) while L2[i]< r do od:

L[i][1]:

end:


#GenM(A,B,p): a random A by B 0-1 matrix  with Pr(1)=p
#(indep.)
GenM:=proc(A,B,p) local a,b:
[seq( [ seq( LD(1-p+p*x,x), b=1..B)], a=1..A)]:
end:








#expands S by one 1
Fatten1:=proc(v,S) local i:

for i from 1 to nops(v) do
 if v[i]=1 and (member(i-1,S)  or member( i+1,S)) and not member(i,S)  then
     RETURN(S union {i}):
 fi:
od:

FAIL:

end:

#Fatten(v,S): inputs a vector of 0-1 v, and a set S of indices 
#with S[i]=1, expands it to the set of all 1's reachable from S
#horiz.
Fatten:=proc(v,S) local S1,S2,S3:

S1:=S:

S2:=Fatten1(v,S1):

while S2<>FAIL and S2<>S1 do
S3:=Fatten1(v,S2):

 S1:=S2:
 S2:=S3:

od:

S1:

end:

  #RO(M): the set of places where M[nops(M)] has
#a 1 reachable from M[1]
RO:=proc(M) local m,S,a,M1,S1:
m:=nops(M):

if m=1 then
S:={}:

for a from 1 to nops(M[1]) do
 if M[1][a]=1 then
    S:=S union {a}:
 fi:
od:
RETURN(S):
fi:

M1:=[op(1..m-1,M)]:

S1:=RO(M1):


S:={}:

for a from 1 to nops(M[m]) do
 if M[m][a]=1  and member(a,S1) then
    S:=S union {a}:
 fi:
od:

Fatten(M[m],S):

end:



#IsGood(M): Given a 0-1 matrix decides whether there is a
#"continuous" path of 1 from the first floor to the nops(M)'s floor
IsGood:=proc(M)
evalb(RO(M)<>{}):
end:

#MC(A,B,K,p):runs IsGood(M) K times applied to random A by B 0-1 matrices
MC:=proc(A,B,p,K) local c,M,i:
c:=0:

for i from  1 to K do
 M:=GenM(A,B,p):
 if IsGood(M) then
    c:=c+1:
 fi:
od:

c/K:

end:

#Wt(M,p): the prob. of the 0-1 matrix M
Wt:=proc(M,p) local i,j,c:

c:=1:

for i from 1 to nops(M) do
 for j from 1 to nops(M[i]) do
   if M[i][j]=1 then
      c:=c*p:
   else
       c:=c*(1-p):
   fi:
 od:
od:

c:

end: