ez:=proc(): print(` M(n,x), KWFun(x,K), MaxEV(n,x), EVKWm(x,m),Per(n), RMR(n,x),L(n), LML(n,x) . CheckKW(n,x) `): end:

with(linalg):

# M(n,x): the 2^(n-1) by 2^(n-1) transfer matrix of H.A. Kramers and G. H. Wannier with x=exp(2*K), the case n=5 is
#Eq. (20) pn p. 260 of their paper "Statistics of the Two-Dimensional Ferromagnet. Part I", Physical Review, v. 60 (Aug. 1, 1941), 252-262
M:=proc(n,x) local i,j,M:
 for i from 0 to 2^(n-1)-1 do
  for j from 0 to 2^(n-1)-1 do
     M[i,j]:=0:
  od:
 od:

for i from 0 to 2^(n-2)-1 do
  M[i,2*i]:=x: 
   M[i,2*i+1]:=1:
od:


for i from 0 to 2^(n-2)-1 do
  M[2^(n-2)+i,2^(n-1)-1-2*i]:=1:
  M[2^(n-2)+i,2^(n-1)-2-2*i]:=1/x:		
od:

[seq([seq(M[i,j],j=0..2^(n-1)-1)], i=0..2^(n-1)-1)]:

end:



#Per(n): the permutation that leads to the matrix Gothic R . Try: Per(5);
Per:=proc(n) local i,gu,T,i1:

for i from 0 to 2^(n-1)-1 do

 gu:=convert(i,base,2):
 gu:=[seq(gu[nops(gu)+1-i1],i1=1..nops(gu))]:
 
 gu:=[0$(n-nops(gu)),op(gu)]:

 gu:=[seq(op([gu[2*i1-1],1-gu[2*i1]]),i1=1..(n-1)/2),gu[nops(gu)]]:


 T[i]:=add(gu[i1]*2^(n-i1),i1=1..nops(gu)):
od:

[seq(T[i]+1,i=0..2^(n-1)-1)]:

end:


#RMR(n,x): the matrix obtained by permuting the rows and columns by Per(n)
RMR:=proc(n,x) local pi,M1,i,j:
pi:=Per(n):

M1:=M(n,x):

[seq([seq(M1[pi[i]][pi[j]],j=1..2^(n-1))],i=1..2^(n-1))]:
end:

#L(n): the matrix L whose n=5 is Eq. (25) in the Kramers-Wannier paper
#paper "Statistics of the Two-Dimensional Ferromagnet. Part I", Physical Review, v. 60 (Aug. 1, 1941), 252-262
L:=proc(n) local gu,mu,i,mu1,i1:
option remember:

if n=1 then
 RETURN([[1]]):
fi:

gu:=[]:
mu:=L(n-1):

for i from 1 to nops(mu) do
 mu1:=mu[i]:
 mu1:=[seq(mu1[i1]/sqrt(2),i1=1..nops(mu1))]:
 gu:=[op(gu),[op(mu1),op(mu1)],[op(mu1),-op(mu1)]]:
od:
gu:
end:


#LML(n,x): the matrix L(n)M(n,x)L[n]. Try: LML(3,x);
LML:=proc(n,x) local gu,mu:
gu:=Matrix(L(n)):
mu:=Matrix(M(n,x)):

multiply(multiply(gu,mu),gu):

end:


#CheckKW(n): checks that LML(n,x) equals (x+1/x)/2 times the transpose of M(n,x*) where x*=(x-1)/(x+1). Try:
#CheckKW(3);
CheckKW:=proc(n) local gu,mu,i,j,x:

gu:=LML(n,x):

gu:=expand([seq([seq(gu[i,j],j=1..2^(n-1))],i=1..2^(n-1))]):

mu:=M(n,x):
mu:=[seq([seq(mu[i][j],i=1..nops(mu))],j=1..nops(mu))]:
mu:=subs(x=(x+1)/(x-1),mu):
mu:=[seq([seq((x-1/x)/2*mu[i][j],j=1..nops(mu[i]))],i=1..nops(mu))]:
mu:=normal(mu):
mu:=expand(mu):

evalb(expand(charpoly(gu,X)-charpoly(mu,X))=0):
end: