Help:=proc(): print(`IsComp(v1,v2), BT1(v) , BT(v) , PerPf(M,N,p),`):
print(` PerPf1(M,N,p,s), BottomStates(M)`): end:

HelpOld:=proc(): print(`PerP(M,N,p), States(M) `): 
print(`IsGoodState(M) `):
end:


HelpAncient:=proc(): print(` LD(p,x) , GenM(A,B,p), Fatten1(V,S) , Fatten(v,S)`): 
print(`RO(M), MC(A,B,p,K) , 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,p,K):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:


###end old stuff
##############start new stuff for April 2, 2009

#PerP(M,N,p): the exact expression (as a poly. in p)
#of the prob. of (simplified) percolation from
#the bottm row to the top row
PerP:=proc(M,N,p) local c,i,v:
c:=0:
for i from 0 to 2^(M*N)-1 do
 v:=convert(i,base,2):
  v:=[op(v),0$(M*N-nops(v))]:
  v:=[seq( [op(N*j+1..N*(j+1),v)],j=0..M-1)]:
 
if IsGood(v) then
  c:=expand(c+Wt(v,p)):
fi:

od:
c:
end:


#States(M): all the vectors of size M
#with {0,1,2} with at least one 2
States:=proc(M) local S,i,v:
option remember:

S:={}:

for i from 0 to 3^M-1 do
v:=convert(i,base,3):
v:=[op(v),0$(M-nops(v))]:

if IsGoodState(v) then

S:=S union {v}:

fi:

od:

S:

end:

IsGoodState:=proc(v) local i: 

if not member(2,convert(v,set)) then
 RETURN(false):
fi:

for i from 1 to nops(v)-1 do
  if {v[i],v[i+1]}={1,2} then
      RETURN(false):
  fi:
od:

true:
end:


###new stuff for Apr. 6, 2009####



#BT1(v): all smallest pair [i,j] such that v[i]=v[i+1]=...=v[j]=2
BT1:=proc(v) local L,i,j:

L:=[]:

for i from 1 to nops(v) while v[i]<>2 do od:

for j from i+1 to nops(v) while v[j]=2 do od:
 

  [i,j-1]:
end:


#BT(v): The list of pairs [i,j] such that v[i]=v[i+1]=...=v[j]=2
BT:=proc(v) local L,i,j,emilie,v1,L1,x:

if v=[] then
 RETURN([]):
fi:

emilie:=BT1(v):

j:=emilie[2]:

if j>nops(v) then
 RETURN([]):
fi:


v1:=[op(j+1..nops(v),v)]:

L1:=BT(v1):

L1:=[seq(L1[x]+[j,j],x=1..nops(L1))]:

[emilie,op(L1)]:

end:


#IsComp(v1,v2) returns true iff state v2 can reside right on top
#of state v1. For example IsComp([0,2],[0,2]); is true
#but IsComp([0,2],[0,1]); is false
IsComp:=proc(v1,v2) local i,L,dan:
 
for i from 1 to nops(v1) do
 if v1[i]=2 and v2[i]=1 then
   RETURN(false):
 fi:
od:


L:=BT(v2):

for dan in L do
 if not member(2,{op(dan[1]..dan[2],v1)}) then
   RETURN(false):
 fi:
od:




true:

end: 

#PerPf(M,N,p): a fast version of PerP(M,N,p) using the transfer matrix
#method, by using States(N) and building up M
PerPf:=proc(M,N,p) local S,s:
S:=States(N):
 
 expand(add(PerPf1(M,N,p,s),s in  S)):

 end:
#Wt2(M,p): the prob. of the 0-1 matrix M
Wt2:=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 or M[i][j]=2 then
      c:=c*p:
   else
       c:=c*(1-p):
   fi:
 od:
od:

c:

end:


#BottomStates(M): all the vectors of size M
#with {0,1,2} with at least one 2
BottomStates:=proc(M) local S,i,v:
option remember:

S:={}:

for i from 0 to 3^M-1 do
v:=convert(i,base,3):
v:=[op(v),0$(M-nops(v))]:

if IsGoodState(v)  and not member(1,v) then

S:=S union {v}:

fi:

od:

S:

end:
















#PerPf1(M,N,p,s)
PerPf1:=proc(M,N,p,s) local c,S,s1,BS:
option remember:

BS:=BottomStates(N):
if M=1 then
  if member(s,BS) then

     RETURN(Wt2([s], p)):
   else
    RETURN(0):
   fi:


fi:

S:=States(N):

c:=0:

for s1 in S do
 if IsComp(s1,s) then
   c:=c+PerPf1(M-1,N,p,s1):
 fi:
od:

expand(c*Wt2([s],p)):

end: