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

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

#Help procedure
Help:=proc():
    print(` MatMulRG(A, B), Strassen(A, B) , RandMat(n, R) `):
end:

##PROBLEM 3##

#MatMultRG(A,B)
#Uses MatMultR(A,B) in order to multiply any two square matrices
#A and B of the same dimension, not necessarily a power of 2. 
MatMultRG:=proc(A, B) local Apad, Bpad, Prod, i, j:
    Apad:=PadZeros(A):
    Bpad:=PadZeros(B):
    Prod:=MatMultR(Apad, Bpad):
    return [seq([seq(Prod[i][j], j=1..nops(A))], i=1..nops(A))]:
end:

##PROBLEM 4##
 
#Strassen(A,B): implements Strassen's algorithm to multiply square
#matrices A and B
Strassen:=proc(A,B)
local Apad, Bpad, Prod, n,i,A11,A12,A21,A22,B11,B12,B21,B22,
        M1,M2,M3,M4,M5,M6,M7,C11,C12,C21,C22:
    if nops(B)<>nops(A) then
        return FAIL:
    fi:
    
    if nops(A)=1 then
        return [[A[1][1]*B[1][1]]]:
    fi:
    
    Apad:=PadZeros(A):
    Bpad:=PadZeros(B):
    
    n:=nops(Apad):

    A11:=[seq([op(1..n/2,Apad[i])],i=1..n/2)]:
    A12:=[seq([op(n/2+1..n,Apad[i])],i=1..n/2)]:
    A21:=[seq([op(1..n/2,Apad[i])],i=n/2+1..n)]:
    A22:=[seq([op(n/2+1..n,Apad[i])],i=n/2+1..n)]:



    B11:=[seq([op(1..n/2,Bpad[i])],i=1..n/2)]:
    B12:=[seq([op(n/2+1..n,Bpad[i])],i=1..n/2)]:
    B21:=[seq([op(1..n/2,Bpad[i])],i=n/2+1..n)]:
    B22:=[seq([op(n/2+1..n,Bpad[i])],i=n/2+1..n)]:



    M1:=Strassen(MatAdd(A11, A22), MatAdd(B11, B22)):
    M2:=Strassen(MatAdd(A21, A22), B11):
    M3:=Strassen(A11, MatAdd(B12, -B22)):#yes, this works...
    M4:=Strassen(A22, MatAdd(B21, -B11)):
    M5:=Strassen(MatAdd(A11, A12), B22):
    M6:=Strassen(MatAdd(A21, -A11), MatAdd(B11, B12)):
    M7:=Strassen(MatAdd(A12, -A22), MatAdd(B21, B22)):
     
    C11:=MatAdd(MatAdd(MatAdd(M1, M4), -M5), M7):
    C12:=MatAdd(M3, M5):
    C21:=MatAdd(M2, M4):
    C22:=MatAdd(MatAdd(MatAdd(M1, -M2), M3), M6):

    Prod:=[seq([op(C11[i]),op(C12[i])],i=1..nops(C11)),
        seq([op(C21[i]),op(C22[i])],i=1..nops(C21))]:
    return [seq([seq(Prod[i][j], j=1..nops(A))], i=1..nops(A))]:
end:

##PROBLEM 5##

#RandMat(n,R): inputs a positive integer n,
#and a large positive integer R,
#and outputs a random n by n matrix with entries
#that are integers between 0 and R
RandMat:=proc(n, R) local r, i, j:
    r:=rand(0..R):
    return [seq([seq(r(), i=1..n)], j=1..n)]:
end:

##PROBLEM 6##
#M:=RandMat(64, 999):
#N:=RandMat(64, 999):
#time(MatMultR(M, N)); returned 14.131
#time(Strassen(M, N)); returned 12.753

#M:=RandMat(128, 9):
#N:=RandMat(128, 9):
#time(MatMultR(M, N)); returned 97.291
#time(Strassen(M, N)); returned 89.673

#There appears to be a small speed-up