#hw4.txt Cole Franks 06/02/1991

###MatMultRG(A,B) calls padzeros to add zeros to A and B so their dimensions are 
#powers of two, uses matmultR to multiply them, and finally removes excess zeros. 


MatMultRG:=proc(A,B) local i,j,C,D,E:

C:=PadZeros(A):

D:=PadZeros(B):

E:=MatMultR(C,D):



[
	seq(
		[
			seq(
				E[i][j],j=1..nops(A)
			)
		], i=1..nops(A)
	)
]:



end:


#####MatMinus

MatMinus:=proc(A):

[
	seq(
		[
			seq(
				-A[i][j], j = 1..nops(A)
			)
		], i = 1..nops(A)
	)
]:


end:

#####Strassen(A,B) implements strassen’s procedure to multiply matrices
#Assuming they are already powers of 2

Strassen:=proc(A,B) local A11,A12, A21, A22, B11, B12, B21, B22, 
C11, C12, C21, C22, M1, M2,
M3, M4, M5, M6, M7, E11, E12, E21, E22, F, n, i, j:

if nops(A) < 2 then

	return(MatMult(A,B)):
fi:

n:=nops(A):

# [ A11, A12]
# [ A21, A22]


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



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




 #   M1 := (A11 + A22) (B11 + B22)
 #   M2 := (A21 + A22) B11
#    M3 := A11 (B12 - B22)
#    M4 := A22 (B21 - B11)
  #  M5 := (A11 + A12) B22
#    M6 := (A21 - A11) (B11 + B12)
#    M7 := (A12 - A22) (B21 + B22) 


M1 := Strassen(MatAdd(A11, A22), MatAdd(B11,B22)):
M2 := Strassen(MatAdd(A21, A22), B11):
M3 := Strassen(A11, MatAdd(B12,MatMinus(B22))):
M4 := Strassen(A22, MatAdd(B21,MatMinus(B11))):
M5 := Strassen(MatAdd(A11, A12), B22):
M6 := Strassen(MatAdd(A21, MatMinus(A11)), MatAdd(B11,B12)):
M7 := Strassen(MatAdd(A12, MatMinus(A22)), MatAdd(B21,B22)):

#add things back together to get the product


C11:= MatAdd(
	M1, MatAdd(
		M4, MatAdd(
			MatMinus(M5), M7
			)
		)
	):

C12:= MatAdd(M3,M5):
C21:=MatAdd(M2,M4):
C22:=MatAdd(
	M1, MatAdd(
		MatMinus(M2), MatAdd(
			M3, M6
			)
		)
	):

#make the returned matrix out of the block matrices

[seq([op(C11[i]),op(C12[i])],i=1..nops(C11)),
seq([op(C21[i]),op(C22[i])],i=1..nops(C21))]:
				


end:


####StrassMult implements Strassen on any square matrices similarly to MatMultRG.

StrassMult:=proc(A,B) local i,j,C,D,E:

C:=PadZeros(A):

D:=PadZeros(B):

E:=Strassen(C,D):



[
	seq(
		[
			seq(
				E[i][j],j=1..nops(A)
			)
		], i=1..nops(A)
	)
]:



end:


###RandMatPos(n, R) outputs a random nXn matrix with random integers between 
### 0 and R

RandMatPos:=proc(n,R) local ra,i,j:
ra:=rand(0..R):
[seq([seq(ra(),j=1..n)],i=1..n)]:
end:

### tests strassen compared to MatMult.

test:=proc(n) local i:
for i from 3 to n do
	print(time(StrassMult(RandMatPos(2^i,1),RandMatPos(2^i,1))),
		time(MatMult(RandMatPos(2^i,1),RandMatPos(2^i,1)))):
od:
end:


### I did not notice any speedup when I ran proc(10). In fact, Strassen did 
#much worse. However, Strassen was comparable to MatMultRG.

#For example, the times for Strassen and MatMult for a 2^9 size 
#square matrix were 3421.225 and 88.17 respectively.