Help:=proc(): print(`PerP(M,N,p), States(M) `): print(`IsGoodState(M) `): end: HelpOld:=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: 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: