#hw4.txt, 03 Feb 2014
#Sowmya Srinivasan


#a few examples from Garvan's booklet pages 91 - end

for i from 1 to 10 do 
   print(i) : 
od : 

with(combinat) : 
choose(10,2); 

with(numtheory) : 
phi(100) ; 


#Read C4.txt for PadZeros and MatMult


#MatMUltRG(A,B) multiplies square matrices A,B of the same dimension
#not necessarily a power of 2
MatMultRG:= proc(A,B)  local ind,APad,BPad,CPad,i,j; 


if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
 print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
 RETURN(FAIL):
fi:


if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
 print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops(A[1])<> nops(B) then
   print (`number of column of`,A,`is not the same as the number of rows of`,B):
   print(`can't multiply these two matrices`):
   RETURN(FAIL):
fi:

ind:= nops(A): 
APad:= PadZeros(A); 
BPad:= PadZeros(B); 
CPad:= MatMultR(APad,BPad):

[seq([seq(CPad[i][j],j=1..ind)],i=1..ind)] :



end:


#MatMinus(A,B) subtracts matrices A,B and is used in Strassen(A,B)
MatMinus:=proc(A,B) local i,j:

if nops(A)<>nops(B) then
 print(`Not the same number of rows`):
 RETURN(FAIL):
fi:

if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
 print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
 print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops(A)>0 and nops(A[1])<>nops(B[1]) then
 print(`Not the same number of columns`):
 RETURN(FAIL):
fi:


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

end:


##Strassen(A,B) uses the Strassen algorithm to multiply matrices A,B
#that are square and have size a power of 2
Strassen:=proc(A,B) local i,n,A11,A12,A21,A22,B11,B12,B21,B22,
M1,M2,M3,M4,M5,M6,M7,C11,C12,C21,C22    : 

n:=nops(A): 

if nops(B)<>n then
 RETURN(FAIL):
fi:

if not type(log[2](n),integer) then
     RETURN(FAIL):
fi:

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

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:=Strassen(MatAdd(A11,A22),MatAdd(B11,B22)): 
M2:=Strassen(MatAdd(A21,A22),B11): 
M3:=Strassen(A11,MatMinus(B12,B22)):
M4:=Strassen(A22,MatMinus(B21,B11)): 

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

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

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

end: 



##StrassenRG(A,B) does the Strassen algorithm for A,B not necessarily of size
#a power of 2
StrassenRG:= proc(A,B)  local ind,APad,BPad,CPad,i,j; 


if nops({seq(nops(A[i]),i=1..nops(A))})<>1 then
 print(`Rows do not all the same lengths, hence`, A, `is not a matrix `):
 RETURN(FAIL):
fi:


if nops({seq(nops(B[i]),i=1..nops(B))})<>1 then
 print(`Rows do not all the same lengths, hence`, B, `is not a matrix `):
 RETURN(FAIL):
fi:

if nops(A[1])<> nops(B) then
   print (`number of column of`,A,`is not the same as the number of rows of`,B):
   print(`can't multiply these two matrices`):
   RETURN(FAIL):
fi:

ind:= nops(A): 
APad:= PadZeros(A); 
BPad:= PadZeros(B); 
CPad:= Strassen(APad,BPad):

[seq([seq(CPad[i][j],j=1..ind)],i=1..ind)] :

end:


#RandMatSq(n,R) outputs a random square matrix of size n 
#with integer entries between 0 and R
RandMatSq:=proc(n,R) local ra,i,j:
ra:=rand(0..R):
[seq([seq(ra(),j=1..n)],i=1..n)]:
end:

#For 64 by 64 matrices the difference between Strassen(A,B) and MatMultR(A,B) 
#does not seem much. For 128 by 128, Strassen(A,B) seems to take much longer
#than MatMultR, though it should be faster.