矩阵快速幂和矩阵的基本操作

参考文章：

http://blog.csdn.net/shiwei408/article/details/8818386

http://www.cnblogs.com/yan-boy/archive/2012/11/29/2795294.html

http://www.cnblogs.com/vongang/archive/2012/04/01/2429015.html

矩阵是线性代数的知识。。。后悔没好好学了。。。

第一部分：矩阵的基础知识

1.结合性 （AB)C=A(BC）．

2.对加法的分配性 （A+B)C=AC+BC，C(A+B）=CA+CB ．

3.对数乘的结合性 k(AB）=（kA)B =A(kB）．

4.关于转置 (AB)'=B'A'．

一个矩阵就是一个二维数组，为了方便声明多个矩阵，我们一般会将矩阵封装一个类或定义一个矩阵的结构体，我采用的是后者。

第二部分：矩阵相乘

若A为n×k矩阵，B为k×m矩阵，则它们的乘积AB(有时记做A·B)将是一个n×m矩阵。

其乘积矩阵AB的第i行第j列的元素为：

举例：A、B均为3*3的矩阵：C=A*B，下面的代码会涉及到两种运算顺序，第一种就是直接一步到位求，第二种就是每次求一列，比如第一次，C00+=a00*b00,C01+=a00*b01……第二次C00+=a00*b10，C01+=a01*b11……以此类推。。。

C00 = a00*b00 + a01*b10 + a02*b20
C01 = a00*b01 + a01*b11 + a02*b21 
C02 = a00*b02 + a01*b12 + a02*b22
C10 = a10*b00 + a11*b10 + a12*b20
C11 = a10*b00 + a11*b11 + a12*b21
C12 = a10*b02 + a11*b12 + a12*b22
C20 = a20*b00 + a21*b10 + a22*b20
C21 = a20*b01 + a21*b11 + a22*b21
C22 = a20*b02 + a21*b12 + a22*b22
C00 = a00*b00 + a01*b10 + a02*b20
C01 = a00*b01 + a01*b11 + a02*b21 
C02 = a00*b02 + a01*b12 + a02*b22

具体该怎么实现两个矩阵相乘呢?

一般会用O(n^3)的方法。。。配合剪枝【添条件，设门槛。。。】，如下：其实主要就是函数 MATRIX mat_multiply (MATRIX a , MATRIX b , int n);

[cpp] view plain copy

1. //O(n^3)算法  
2. #include <iostream>  
3. #include <cstdio>  
4. #include <cstring>  
5. #include <cstdlib>  
6. #include <cmath>  
7. #include <algorithm>  
8. using namespace std;  
9. #define LL __int64  
10. #define MOD 1000  
11. typedef struct MATRIX  
12. {  
13.     int mat[50][50];  
14. }MATRIX;  
15.   
16. MATRIX mat_multiply (MATRIX a,MATRIX b,int n)  
17. {  
18.     MATRIX c;   //c[i][j]= Σ a[i][k]*b[k][j]  
19.     memset(c.mat,0,sizeof(c.mat));      
20. /* 
21.     for(int i=0;i<n;i++)    //a矩阵一行一行往下 
22.         for(int j=0;j<n;j++)    //b矩阵一列一列往右 
23.             for(int k=0;k<n;k++)    //使a矩阵 第i行第k个元素 乘以 b矩阵 第j列第k个元素 
24.                 if(a[i][k] && b[k][j])    //剪枝（添条件，设门槛），提高效率，有一个是0，相乘肯定是0 
25.                     c.mat[i][j]+=a.mat[i][k]*b.mat[k][j]; 
26. */  
27.   
28. //上面也是可以的，但是下面的剪枝更好一些，效率更高一些，但是运算顺序有点难想通，，，  
29. //上面就是C[i][j]一次就求出来，下面就是每次c[i][j]求出一项【就是上面红体字，每次各求一列】  
30.   
31.     for(int k=0;k<n;k++)    //这个可以写到前面来，  
32.         for(int i=0;i<n;i++)  
33.             if(a.mat[i][k])     //剪枝：如果a.mat[i][k]是0，就不执行了  
34.                 for(int j=0;j<n;j++)  
35.                     if(b.mat[k][j])     //剪枝：如果b.mat[i][k]是0，就不执行了  
36.                     {  
37.                         c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];  
38.                         if(c.mat[i][j]>=MOD)    //这个看实际需求，要不要取模  
39.                             c.mat[i][j]%=MOD;   //取模的复杂度比较高，所以尽量减少去模运算，添加条件，只有当大于等于MOD的时候才取余  
40.                     }  
41.     return c;  
42. }  
43.   
44. int main()<span style="white-space:pre">  </span>//这个只是用来测试用的。。。  
45. {  
46.     int n;  
47.     MATRIX A,B,C;  
48.   
49.     memset(A.mat,0,sizeof(A.mat));  
50.     memset(B.mat,0,sizeof(B.mat));  
51.     memset(C.mat,0,sizeof(C.mat));  
52.   
53.     scanf("%d",&n);     //矩阵规模，这里是方阵，行数等于列数  
54.   
55.     for(int i=0;i<n;i++)    //初始化A矩阵  
56.         for(int j=0;j<n;j++)  
57.             scanf("%d",&A.mat[i][j]);  
58.   
59.     for(int i=0;i<n;i++)    //初始化B矩阵  
60.         for(int j=0;j<n;j++)  
61.             scanf("%d",&B.mat[i][j]);  
62.   
63.     C=mat_multiply (A,B,n);  
64.   
65.     for(int i=0;i<n;i++)    //打印C矩阵  
66.     {  
67.         for(int j=0;j<n;j++)  
68.             printf("%3d",C.mat[i][j]);  
69.         printf("\n");  
70.     }  
71.     return 0;  
72. }  

第三部分：矩阵快速幂

神马是幂？【很多时候会被高大上的名字吓到。。。导致学习效率降低。。。其实没辣么可怕，很简单！！！】

幂又称乘方。表示一个数字乘若干次的形式,如n个a相乘的幂为a^n ，或称a^n为a的n次幂。a称为幂的底数，n称为幂的指数。——引自.度娘百科

这类题，指数都是很大很大很大很大很大很大很大的。。。霸王硬上弓的话，很容易超时的 T_T 。。。所以得快速幂→_→

学过之后发现，其实矩阵快速幂 的核心思想跟 以前学过的快速幂取模非常非常相似，只是矩阵乘法需要另外写个函数，就是上面那个代码。。。

【一会去写快速幂取模的专题，棒！】

快速幂的思路就是：

设A为矩阵，求A的N次方，N很大，1000000左右吧。。。

先看小一点的，A的9次方

A^9

= A*A*A*A*A*A*A*A*A  【一个一个乘，要乘9次】

= A*(A*A)*(A*A)*(A*A)*(A*A)【保持格式的上下统一，所以加上这句】

 = A*(A^2)^4 【A平方后，再四次方，还要乘上剩下的一个A，要乘6次】

= A*((A^2)^2)^2【A平方后，再平方，再平方，还要乘上剩下的一个A，要乘4次】

也算是一种二分思想的应用吧，1000000次幂，暴力要乘1000000次，快速幂就只要（log2底1000000的对数） 次，大约20次。。。这。。。我没错吧。。。

但是因为是矩阵，矩阵乘法的复杂度是O(n^3)。。。所以快速幂的复杂度是O(n^3 * logn)

上代码！矩阵乘法的代码和上面一样，添加了mat_quick_index(MATRIX a,int N,int n)函数，主函数做了些许修改，以便检验。。。

[cpp] view plain copy

1. //矩阵快速幂  
2. #include <iostream>  
3. #include <cstdio>  
4. #include <cstring>  
5. #include <cstdlib>  
6. #include <cmath>  
7. #include <algorithm>  
8. using namespace std;  
9. #define LL __int64  
10. #define MOD 1000  
11. typedef struct MATRIX  
12. {  
13.     int mat[50][50];  
14. }MATRIX;  
15.   
16. MATRIX mat_multiply (MATRIX a,MATRIX b,int n)  
17. {  
18.     MATRIX c;   //c[i][j]= Σ a[i][k]*b[k][j]  
19.     memset(c.mat,0,sizeof(c.mat));  
20. /* 
21.     for(int i=0;i<n;i++)    //a矩阵一行一行往下 
22.         for(int j=0;j<n;j++)    //b矩阵一列一列往右 
23.             for(int k=0;k<n;k++)    //使a矩阵 第i行第k个元素 乘以 b矩阵 第j列第k个元素 
24.                 if(a[i][k] && b[k][j])    //剪枝（添条件，设门槛），提高效率，有一个是0，相乘肯定是0 
25.                     c.mat[i][j]+=a.mat[i][k]*b.mat[k][j]; 
26. */  
27.   
28. //上面也是可以的，但是下面的剪枝更好一些，效率更高一些，但是运算顺序有点难想通，，，  
29. //上面就是C[i][j]一次就求出来，下面就是每次c[i][j]求出一项【就是上面红体字，每次各求一列】  
30.   
31.     for(int k=0;k<n;k++)    //这个可以写到前面来，  
32.         for(int i=0;i<n;i++)  
33.             if(a.mat[i][k])     //剪枝：如果a.mat[i][k]是0，就不执行了  
34.                 for(int j=0;j<n;j++)  
35.                     if(b.mat[k][j])     //剪枝：如果b.mat[i][k]是0，就不执行了  
36.                     {  
37.                         c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];  
38.                         if(c.mat[i][j]>=MOD)    //这个看实际需求，要不要取模  
39.                             c.mat[i][j]%=MOD;   //取模的复杂度比较高，所以尽量减少去模运算，添加条件，只有当大于等于MOD的时候才取余  
40.                     }  
41.     return c;  
42. }  
43.   
44. MATRIX mat_quick_index(MATRIX a,int N,int n)  
45. {  
46.     MATRIX E;   //单位矩阵，就像数值快速幂里，把代表乘积的变量初始化为1  
47.   
48. //    memset(E.mat,0,sizeof(E.mat));  //置零，单位矩阵除了主对角线都是1，其他都是0  
49. //    for(int i=0;i<n;i++)    //初始化单位矩阵【就是主对角线全是1】  
50. //            E.mat[i][i]=1;  
51.   
52.     for(int i=0;i<n;i++)  
53.         for(int j=0;j<n;j++)  
54.             E.mat[i][j]=(i==j); //酷炫*炸天的初始化！！！  
55.   
56.     while(N>0)  
57.     {  
58.         if(N & 1)  
59.             E=mat_multiply(E,a,n);  
60.         N>>=1;  
61.         a=mat_multiply(a,a,n);  
62.     }  
63.     return E;  
64. }  
65.   
66. int main()  
67. {  
68.     int n,N;    //n为矩阵（方阵）规模，几行，N为指数  
69.     MATRIX A,C;  
70.   
71.     memset(A.mat,0,sizeof(A.mat));  
72.     memset(C.mat,0,sizeof(C.mat));  
73.   
74.     scanf("%d",&n);     //矩阵规模，这里是方阵，行数等于列数  
75.   
76.     for(int i=0;i<n;i++)    //初始化A矩阵  
77.         for(int j=0;j<n;j++)  
78.             scanf("%d",&A.mat[i][j]);  
79.   
80.     scanf("%d",&N);  
81.     C=mat_quick_index(A,N,n);  
82.   
83.     for(int i=0;i<n;i++)    //打印C矩阵  
84.     {  
85.         for(int j=0;j<n;j++)  
86.             printf("%3d",C.mat[i][j]);  
87.         printf("\n");  
88.     }  
89.     return 0;  
90. }  

矩阵快速幂暂时就是这么多啦~~~【如果有新的东西，会更新的~~~】

主要就是明白它的原理，然后根据实际情况进行修改代码即可~~~！！！

最后 推荐几道矩阵快速幂的题目吧：

POJ：3070、3150、3233、3735

好啦~~~写了好久的感觉，但是觉得自己对矩阵快速幂也有了更加深入的理解了~~~！！！

转自https://blog.csdn.net/u013795055/article/details/38599321