用Strassen算法实现矩阵乘法

#coding=utf-8

'''第10题：用Strassen的思想实现两个n×n矩阵的乘法'''

import math
from numpy import *

#check函数用来检查num是否是2的整数次幂
def check(num):
i = 1;
while (True):
if (i > num):
return False
if (i == num):
return True
i = i * 2

#Multiply函数即为实现方法
def Matrix_Multiply(A,B):
#由于本题只考虑两个n×n矩阵的乘法，且n为2的整数次幂的情况，故做一些判断去掉不合法的输入
if A.shape!=B.shape:
print '本题只考虑两个n×n矩阵的乘法，且n为2的整数次幂的情况'
return None
if A.shape[0]!=A.shape[1]:
print '本题只考虑两个n×n矩阵的乘法，且n为2的整数次幂的情况'
return None
if check(A.shape[0])==False:
print '本题只考虑两个n×n矩阵的乘法，且n为2的整数次幂的情况'
return None
n = A.shape[0]
result = zeros([n,n])

if(n == 2):
#最基本的情况
A11=A[0,0]
A12=A[0,1]
A21=A[1,0]
A22=A[1,1]
B11=B[0,0]
B12=B[0,1]
B21=B[1,0]
B22=B[1,1]

P1 = A11*(B12-B22)
P2 = (A11+A12)*B22
P3 = (A21+A22)*B11
P4 = A22*(B21-B11)
P5 = (A11+A22)*(B11+B22)
P6 = (A12-A22)*(B21+B22)
P7 = (A11-A21)*(B11+B12)
else:
#输入复杂的情况就采取divide-and-conquer策略
A11=A[:n/2,:n/2]
A12=A[:n/2,n/2:]
A21=A[n/2:,:n/2]
A22=A[n/2:,n/2:]

B11=B[:n/2,:n/2]
B12=B[:n/2,n/2:]
B21=B[n/2:,:n/2]
B22=B[n/2:,n/2:]
#combine
P1 = Matrix_Multiply(A11,B12-B22)
P2 = Matrix_Multiply(A11+A12,B22)
P3 = Matrix_Multiply(A21+A22,B11)
P4 = Matrix_Multiply(A22,B21-B11)
P5 = Matrix_Multiply(A11+A22,B11+B22)
P6 = Matrix_Multiply(A12-A22,B21+B22)
P7 = Matrix_Multiply(A11-A21,B11+B12)

result[:n/2,:n/2] = P4 + P5 + P6 - P2 #C11
result[:n/2,n/2:] = P1 + P2 #C12
result[n/2:,:n/2] = P3 + P4 #C21
result[n/2:,n/2:] = P1 + P5 - P3 - P7 #C22

return result

'''demo'''
print '矩阵A='
a = array([[1.1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
print a
print '计算A×A'
print ''
result = Matrix_Multiply(a,a)
print '本题所实现的方法结果是：'
print result
print '正确结果应为：'
print dot(a,a)