#Nathan Fox
#Homework 28
#I give permission for this work to be posted online

#Read procedures from class/last homework
read(`hw27.txt`):
read(`C28.txt`):

#Help procedure
Help:=proc():
    print(` Renormalize(M) , DrawRenormalization(k, z, K, h) `):
    print(` SMagS1(DI, MAT) , SMagS(DI, K1) `):
end:

##PROBLEM 2##

#Renormalize(M): inputs a square matrix (given as a list of lists)
#with 3^k rows and 3^k columns (it should return FAIL if it is not
#of such a size) and outputs the 3^(k-1) by 3^(k-1) matrix obtained
#by splitting it into 3^(k-1) by 3^(k-1) blocks of 3 by 3 matrices,
#and replacing each such block by sign of the sum of the entries.
Renormalize:=proc(M) local ret, n, i, j, k, l:
    if nops(M) <> nops(M[1]) or nops(M) < 3 or not type(log[3](nops(M)), integer) then
        return FAIL:
    fi:
    
    n:=nops(M)/3:
    ret:=[[0$n]$n]:
    for i from 1 to n do
        for j from 1 to n do
            ret[i][j]:=sign(add(add(M[3*i+k][3*j+l], l=-2..0), k=-2..0)):
        od:
    od:
    return ret:
end:

##PROBLEM 3##

#Draw renormalization matrices starting from 3^k by 3^k down to 1
DrawRenormalization:=proc(k, z, K, h) local M, ret, hh:
    M:=Ising(3^k, 3^k, z, K):
    ret:=NULL:
    hh:=h:
    while true do
        ret:=ret, display(DrawMat(M, hh)
        #, textplot([0, -2*h*3^k, sprintf("size=%d", nops(M))])):
        )$5:
        if nops(M) = 1 then
            break:
        fi:
        M:=Renormalize(M):
        hh:=3*h:
    od:
    return display(ret, insequence=true, axes=none):
end:

#DrawRenormalization(4, 1/2, 5000, 3); can be found at hw27_half.gif
#DrawRenormalization(4, 100, 5000, 3); can be found at hw27_100.gif
#These pretty much reproduce Wilson's results

##PROBLEM 4##

#There are too many to find in any reasonable amount of time
#The answer is at least 6710, as that's how many there are using
#only the first 300 words

##PROBLEM 5##

#Symmetric version
SMagS1:=proc(DI, MAT) local S1, I1, J1, W, MAT1, Hope:
    S1:={}:
    Hope:=BEG(DI, nops(MAT) + 1):
    for I1 from 1 to nops(DI) do
        W:=DI[I1]:
        if nops(W) = nops(MAT[1]) and {seq(evalb(W[J1] = MAT[J1][nops(MAT)+1]), J1=1..nops(MAT))} = {true} then
            MAT1:=[op(MAT), W]:
            if {seq([seq(MAT1[I1][J1], I1=1..nops(MAT1))], J1=1..nops(MAT1[1]))} subset Hope then
                S1:=S1 union {MAT1}:
            fi:
        fi:
    od:
    return S1:
end:

#SMagS(DI, K1): outputs symmetric matrices of letters, where the
#rows are identical to the corresponding columns.
SMagS:=proc(DI, K1) local S, I1, MAT:
    S:={seq([DI[I1]] , I1=1..nops(DI))}:
    for I1 from 2 to K1 do
        S:={seq(op(SMagS1(DI, MAT)), MAT in S)}:
    od:
    return S:
end:

#This one actually finished on the whole of E3.  There are 58872
#such squares