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

#Read procedures from class
read(`C5.txt`):

#List Tools
with(ListTools):

#Help procedure
Help:=proc():
    print(` RPower(A, D1) , Onsm(m,z) , NumerOns(m,reso) `):
end:

##PROBLEM 1##

#RPower(A, D1)
#
#INPUT:
#A: a square matrix (given as a list of lists)
#D1: a positive integer
#
#OUTPUT:
#An approximation for the dominant eignevalue with error that is
#believed to be less than 10^-D1, by do-loping over k and
#quiting when the difference between RPow(A,k,x0) and
#RPow(A,k+10,x0) is less than 10^-(D1+2)
#x0 starts as all 1s (everything else was too slow)
#
#NOTE:
#This implementation is stupid.  A better search algorithm should
#be used, but we are required to use RPow(A,k,x0) as a subroutine
#This implementation is also stupid since it may never terminate
#if A doesn't have a dominant eigenvalue
RPower:=proc(A, D1) local n, x0, k, rp, last:
    n:=nops(A):
    x0:=[[1]$n]:
    k:=10:
    last:=infinity:
    while true do
        rp:=RPow(A, k, x0):
        if abs(rp - last) < 10^(-D1-2) then
            return rp:
        fi:
        last:=rp:
        k:=k+10:
    od:
end:

##PROBLEM 2##

#Onsm(m,z)
#
#INPUT:
#a positive integer m
#a number (or variable) z
#
#OUTPUT:
#the famous 2^m by 2^m Onsager matrix defined as follows.
#
#Both rows and columns are labeled by vectors of length m whose
#entries are {-1,1}, where the i-th row (and column) corresponds
#to the {-1,1} vector obtained by converting the integer i to
#binary, and padding with 0's in front to make it of length m,
#and then replacing 0 by -1. Then
#(if u and v are vectors as above) 
#M[u,v]:=z^(u[1]*v[1]+u[2]*v[2]+ ... + u[m]*v[m])*
#z^(u[1]*u[2]+u[2]*u[3]+ ... + u[m-1]*u[m]+u[m]*u[1])
Onsm:=proc(m,z) local vectorize, circulate, i, j:
    vectorize:=proc(n) local L, k:
    option remember:
        L:=Reverse(convert(n, base, 2)):
        L:=[0$(m-nops(L)), op(L)]:
        L:=[seq(2*L[k]-1, k=1..nops(L))]:
        return [L]:
    end:
    circulate:=proc(v)
        return [[op(2..,v[1]), v[1][1]]]:
    end:
    return [seq([seq(z^(MatMult(vectorize(i),
        Tra(vectorize(j)))[1][1])*z^(MatMult(vectorize(i),
        Tra(circulate(vectorize(i))))[1][1]), j=0..2^m-1)],
        i=0..2^m-1)]:
end:

##PROBLEM 3##

#NumerOns(m,reso)
#
#INPUT:
#a positive integer m
#a floating point value reso
#
#OUTPUT:
#the numerical values of fm(z) for z between 0.2 and 3 in
#intervals of reso
#
#NOTE:
#fm(z)=log of largest eigenvalue of Osnm(m,z) divided by m
NumerOns:=proc(m, reso) local z, ret, D1:
    ret:=[]:
    D1:=5:
    for z from 0.2 by reso to 3 do
        ret:=[op(ret), [z, log(RPower(evalf(Onsm(m, z)), D1))/m]]:
    od:
    return ret:
end:

#fm has a minimum around z=1 for all m I tested