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: