矩阵乘法算法

转载自：http://blog.csdn.net/q3498233/article/details/5786180

矩阵运算是属于线性代数里的一个重要内容，上学期学完后只觉得矩阵能解线性方程，不过高中的时候听说过矩阵能优化常系数递推以及将坐标上的点作线性变换，于是找了些资料研究了一下，并把许多经典题以及HDU shǎ崽大牛 总结的矩阵乘法的题目[1] 、[2] 和开设的矩阵乘法DIY Contest 给做完了，感觉收获颇丰。
    一个矩阵就是一个二维数组，为了方便声明多个矩阵，我们一般会将矩阵封装一个类或定义一个矩阵的结构体，我采用的是后者：

[c-sharp]  view plain copy print ?

1. struct Mat  
2. {  
3.     int mat[Max][Max];  
4. }  

 

    最特殊的矩阵应该就是单位矩阵E了，它的对角线的元素为1，非对角线元素为0。

 

一般矩阵乘法采用朴素的O(n^3)的算法：

 

[c-sharp]  view plain copy print ?

1. Mat operator*(Mat a,Mat b)    
2. {    
3.     int i,j,k;    
4.     Mat c;    
5.     for (i=0;i<len;i++)    
6.     {    
7.         for (j=0;j<len;j++)    
8.         {    
9.             c.mat[i][j] = 0;    
10.             for(k=0;k<len;k++)    
11.                 c.mat[i][j]+=(a.mat[i][k]*b.mat[k][j])%MOD;    
12.         }    
13.     }    
14.     return c;    
15. }   

 

 

矩阵加法就是简单地将对应的位置的两个矩阵的元素相加：

[c-sharp]  view plain copy print ?

1. Mat operator+(Mat a,Mat b)    
2. {    
3.     Mat c;    
4.     int i,j;    
5.     for (i=0;i<len;i++)    
6.     {    
7.         for(j=0;j<len;j++)    
8.             c.mat[i][j] = (a.mat[i][j]+b.mat[i][j])%MOD;    
9.     }    
10.     return c;    
11. }   

 

    在ACM的题目中，我们一般考虑的是n阶方阵之间的乘法以及n阶方阵与n维向量（把向量看成n×1的矩阵）的乘法。矩阵乘法最重要的性质就是满足结合律，同时它另一个很重要的性质就是不 满足交换率，这保证了矩阵的幂运算满足快速幂取模 （A^x % MOD）算法：

假设k = 27，则k的二进制表示为11011，所以 

，可以看出：k的二进制的每一位矩阵A都要平方，在k二进制为1的位：末矩阵×平方后的A，在k二进制为0的位则末矩阵×E（单位矩阵），即不变。代码如下：

[c-sharp]  view plain copy print ?

1. Mat operator^(Mat a,int x)    
2.  {    
3.      Mat p = e,q = a;    
4.      while (x)    
5.      {    
6.          if(x&1)    
7.              p = p*q;    
8.          x>>=1;    
9.          q = q*q;    
10.      }    
11.      return p;    
12.  }  

 

许多题目还要求S = A + A ² + A ³ + … + A^k .。其实再作一次二分即可：只需计算log(n)个A的幂即可。

[c-sharp]  view plain copy print ?

1. Mat solve(Mat a,int p)    
2.  {    
3.      if(p==1)    
4.          return a;   
5.      /*如果p为奇数，则对p-1进行二分，a^p+二分结果*/   
6.      /* 
7.           A^1+A^2+A^3+A^4+A^5+A^6+A^7 = (A^1+A^2+A^3)+ A^3*(A^1+A^2+A^3)+A^7=A^7+solve(A,p-1) 
8.      */  
9.      else if(p&1)    
10.          return (a^p)+solve(a,p-1);    
11.      /*p为偶数，则直接二分*/  
12.      /* 
13.      A^1+A^2+A^3+A^4+A^5+A^6 = (A^1+A^2+A^3)+ A^3*(A^1+A^2+A^3) = (A^3+e)*(A^1+A^2+A^3) 
14.      推广即可得到当p为偶数时   A^1+A^2+A^3+……+A^p = (A^(p/2)+e) * solve(A,p/2); 
15.      */  
16.      else    
17.          return ((a^(p>>1))+e)*solve(a,p>>1);    
18.  }    

 

    矩阵乘法算法基本就是这样，欢迎大家指出更好的方法或者提出不足之处