Help:=proc(): print(`En(A), PF(m,n,z) , RealAve1(m,n,z,w,p,q)`): 
print(`RealAve2(m,n,z,w,p1,q1,p2,q2), Cor(m,n,z,w,p1,q1,p2,q2) `):
end:

#En(A): the sum of the products of all nearest-neighbors of 
#the matrix A given as a list of lists
En:=proc(A) local i,j:

 add(add(A[i][j]*A[i][j+1],j=1..nops(A[i])-1),i=1..nops(A))+

add(add(A[i-1][j]*A[i][j],i=2..nops(A)),j=1..nops(A[1])):

end:

#PF(m,n,z)
PF:=proc(m,n,z) local s,c,ab:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab):
od:

s:

end:


#Ex: evaluate f(z)=lim (m,n->infinity) PF(m,n,z)^(1/(m*n))
#(equivalenty: limit of log(PF(m,n,z))/(m*n) )


#PF(m,n,z)
gPF:=proc(m,n,z,w) local s,c,ab,ab1:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab1)*w^convert(ab,`+`):
od:

s:

end:

#Ex: evaluate g(z,w)=lim (m,n->infinity) PF(m,n,z,w)^(1/(m*n))
#(equivalenty: limit of log(PF(m,n,z,w))/(m*n) )


#RealAve1(m,n,z,w,p,q): the (weighted) average of A[i][j]
#amongsts all m by n matrices with {-1,1} with Boltzmann weight
RealAve1:=proc(m,n,z,w,p,q) local t,s,ab,ab1,i,j,c:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab1)*w^convert(ab,`+`)*ab1[p][q]:
od:

s/gPF(m,n,z,w):

end:




#RealAve2(m,n,z,w,p1,q1,p2,q2): the (weighted) average of A[p1][q1]*A[p2][q2]
#amongsts all m by n matrices with {-1,1} with Boltzmann weight
RealAve2:=proc(m,n,z,w,p1,q1,p2,q2) local t,s,ab,ab1,i,j,c:

s:=0:

for c from 0 to 2^(m*n)-1 do
 
 ab:=convert(c,base,2):
 ab:=[op(ab),0$(m*n-nops(ab))]:
 ab:=subs(0=-1,ab):
 ab1:=[seq([op((i-1)*n+1..i*n,ab)],i=1..m)]:
  
 s:=s+z^En(ab1)*w^convert(ab,`+`)*ab1[p1][q1]*ab1[p2][q2]:
od:

s/gPF(m,n,z,w):

end:

#E(XY)-E(X)E(Y)=E((X-E(X))(Y-E(Y))
Cor:=proc(m,n,z,w,p1,q1,p2,q2):

RealAve2(m,n,z,w,p1,q1,p2,q2)-
RealAve1(m,n,z,w,p1,q1)*RealAve1(m,n,z,w,p2,q2):

end: