#Code by Simao M. HERDADE, March 10, 2011

Help:=proc(): print( ` CRMP(A,B)  `): end:



##SubM(A,Li,Lj) gives the submatrix of A forme by the rows  
## of A indexed by Li and the collumns of A indexed by Lj
##
SubM:=proc(A,Li,Lj) local L,i,j,Laux:
L:=[]:
for i from 1 to nops(Li) do 
 	Laux:=[]:
	for j from 1 to nops(Lj) do
		Laux:=[op(Laux),A[Li[i]][Lj[j]]]:
	od:
	L:=[op(L),Laux]:
od:
end:



##FromBlocks(A,B,C,E) creates a square matrix from the
## 4 square matrices A,B,C and E of the same size
##
FromBlocks:=proc(A,B,C,E) local n,M,Aux,i1,j1,j2,i2,j3,j4:
n:=nops(A):
M:=[]:

for i1 from 1 to n do:
	Aux:=[]:
	for j1 from 1 to n do
	Aux:=[op(Aux),A[i1][j1]]:
	od:
	for j2 from 1 to n do
	Aux:=[op(Aux),B[i1][j2]]:
	od:
	M:=[op(M),Aux]:
od:

for i2 from 1 to n do:
	Aux:=[]:
	for j3 from 1 to n do
	Aux:=[op(Aux),C[i2][j3]]:
	od:
	for j4 from 1 to n do
	Aux:=[op(Aux),E[i2][j4]]:
	od:
	M:=[op(M),Aux]:
od:
M:
end:




##Aug(A,g) ads 0's to rows and collums of A to make it
## a square matrix of size g
##
Augm:=proc(A,g) local n,m,Z,M,i,j,l:
n:=nops(A):
m:=nops(A[1]):
Z:=[seq(0,i=1..g)]:
M:=A:

for i from 1 to n do
	M[i]:=[op(M[i]),seq(0,j=m+1..g)]:
	
od:

for l from nops(A)+1 to g do
	M:=[op(M),Z]:
	
od:
M;
end:



##Mad(A,B) does the Matrix addition A+B
##
Mad:=proc(A,B) local n,L,i:
n:=nops(A):
L:=[]: 
for i from 1 to n do
 L:=[op(L),A[i]+B[i]]:
od:
L:
end:


##Dif(A,B) does the Difference A-B
##
Dif:=proc(A,B) local n,L,i:
n:=nops(A):
L:=[]: 
for i from 1 to n do
 L:=[op(L),A[i]-B[i]]:
od:
L:
end:


##CRMP(A,B): RECURSIVE clever multiplication of A and B
##where A and B are arbitrary matrices, using STRASSEN ALGORITHM
##
CRMP:=proc(A,B) local g,k,C,E,s1,s2,i,j, A1,A2,A3,A4,B1,B2,B3,B4,M1,M2,M3,M4,M5,M6,M7,C1,C2,C3,C4,M,o1,o2:


g:=max(nops(A),nops(A[1]),nops(B),nops(B[1])):

if g=1 then
RETURN([[A[1][1]*B[1][1]]]):
fi:


k:=ceil(log[2](g)):

if g< 2^k then
C:=Augm(A,2^k):
E:=Augm(B,2^k):
else
C:=A:
E:=B:
fi:


s1:={seq(i,i=1..2^(k-1))}:
s2:={seq(j,j=2^(k-1)+1..2^k)}:


A1:=SubM(C,s1,s1):
A2:=SubM(C,s1,s2):
A3:=SubM(C,s2,s1):
A4:=SubM(C,s2,s2):

B1:=SubM(E,s1,s1):
B2:=SubM(E,s1,s2):
B3:=SubM(E,s2,s1):
B4:=SubM(E,s2,s2):


##M1=(A1+A4)*(B1+B4), M2=(A3+A4)*B1, M3=A1*(B2-B4), M4=A4*(B3-B1),
##M5=(A1+A2)*B4, M6=(A3-A1)*(B1+B2), M7=(A2-A4)*(B3+B4)
##
M1:=CRMP(Mad(A1,A4),Mad(B1,B4)): 
M2:=CRMP(Mad(A3,A4),B1):
M3:=CRMP(A1,Dif(B2,B4)):
M4:=CRMP(A4,Dif(B3,B1)):
M5:=CRMP(Mad(A1,A2),B4):
M6:=CRMP(Dif(A3,A1),Mad(B1,B2)): 
M7:=CRMP(Dif(A2,A4),Mad(B3,B4)): 


##construct blocks C1,C2,C3,C4
C1:=Dif(Mad(Mad(M1,M4),M7),M5):
C2:=Mad(M3,M5):
C3:=Mad(M2,M4):
C4:=Dif(Mad(Mad(M1,M3),M6),M2):


##construct M from Blocks
M:=FromBlocks(C1,C2,C3,C4):


#get rid of all the extra 0's
o1:={seq(i,i=1..nops(A))}:
o2:={seq(j,j=1..nops(B[1]))}:
SubM(M,o1,o2):

end: