矩阵快速幂初探------由浅入深
矩阵快速幂基本原理：http://www.cnblogs.com/yan-boy/archive/2012/11/29/2795294.html
基础应用：快速Fabonacci算法。
题目链接：http://acm.nyist.net/JudgeOnline/problem.php?pid=148
只要看懂矩阵快速幂的基本原理和这个题目的意思，基本上就可以做了。
核心问题就是构造基本矩阵，设Fabonacci数列为F[n]，需要构造的矩阵为A：
好吧，为了不造成误导，此题中A和Basic我还是用两个矩阵来表示，大家写的时候可以合并掉，这里也贴一下合并后的代码。
1 #include<stdio.h> 2 #include<iostream> 3 using namespace std; 4 #include<queue> 5 #include<math.h> 6 #include<algorithm> 7 #include<string.h> 8 #include<stdlib.h> 9 #include<vector> 10 #include<set> 11 #include<map> 12 13 #define repA(p,q,i) for( int (i)=(p); (i)!=(q); ++(i) ) 14 #define repAE(p,q,i) for( int (i)=(p); (i)<=(q); ++(i) ) 15 #define repD(p,q,i) for( int (i)=(p); (i)!=(q); --(i) ) 16 #define repDE(p,q,i) for( int (i)=(p); (i)>=(q); --(i) ) 17 #define range 1010 18 19 20 struct Matrix{ 21 int v[3][3] ; 22 }; 23 24 Matrix basic; 25 Matrix mtPow( int n ); 26 Matrix mtMultiply( Matrix m, Matrix n ); 27 28 int main() 29 { 30 basic.v[1][1]=0; basic.v[1][2]=1; 31 basic.v[2][1]=1; basic.v[2][2]=1; 32 33 int n; 34 while( scanf("%d",&n) != EOF ) 35 { 36 if(n == -1) break; 37 if(n==0) printf("0\n"); 38 else 39 { 40 Matrix result = mtPow(n) ; 41 printf("%d\n", result.v[1][2] ); 42 } 43 } 44 return 0; 45 } 46 47 Matrix mtPow( int n ) 48 { 49 Matrix result,temp; 50 if(n == 1) return basic ; 51 else 52 { 53 temp = mtPow(n/2) ; 54 result = mtMultiply( temp, temp ) ; 55 } 56 if(n % 2 != 0) 57 return mtMultiply(result , basic) ; 58 else return result ; 59 } 60 61 Matrix mtMultiply( Matrix m, Matrix n ) 62 { 63 Matrix result ; 64 result.v[1][1] = (m.v[1][1]*n.v[1][1] + m.v[1][2]*n.v[2][1]) % 10000 ; 65 result.v[1][2] = (m.v[1][1]*n.v[1][2] + m.v[1][2]*n.v[2][2]) % 10000 ; 66 result.v[2][1] = (m.v[2][1]*n.v[1][1] + m.v[2][2]*n.v[2][1]) % 10000 ; 67 result.v[2][2] = (m.v[2][1]*n.v[1][2] + m.v[2][2]*n.v[2][2]) % 10000 ; 68 return result ; 69 }
代码虽然长了一点，但不至于造成误导，让大家误以为Basic就是A，其实两者完全是两个含义。
Basic是DP的初始状态，而A代表的其实是DP的递推式。
1 #include<stdio.h> 2 #include<iostream> 3 using namespace std; 4 #include<queue> 5 #include<math.h> 6 #include<algorithm> 7 #include<string.h> 8 #include<stdlib.h> 9 #include<vector> 10 #include<set> 11 #include<map> 12 13 #define repA(p,q,i) for( int (i)=(p); (i)!=(q); ++(i) ) 14 #define repAE(p,q,i) for( int (i)=(p); (i)<=(q); ++(i) ) 15 #define repD(p,q,i) for( int (i)=(p); (i)!=(q); --(i) ) 16 #define repDE(p,q,i) for( int (i)=(p); (i)>=(q); --(i) ) 17 #define range 1010 18 19 20 struct Matrix{ 21 int v[3][3] ; 22 }; 23 24 Matrix basic,A; 25 Matrix mtPow( int n ); 26 Matrix mtMultiply( Matrix m, Matrix n ); 27 28 int main() 29 { 30 basic.v[1][1]=0; basic.v[1][2]=1; 31 basic.v[2][1]=1; basic.v[2][2]=1; 32 A.v[1][1]=0; A.v[1][2]=1; 33 A.v[2][1]=1; A.v[2][2]=1; 34 int n; 35 while( scanf("%d",&n) != EOF ) 36 { 37 if(n == -1) break; 38 if(n==0) printf("0\n"); 39 else if(n==1) printf("1\n"); 40 else 41 { 42 Matrix result = mtMultiply( mtPow(n-1) , basic ); 43 printf("%d\n", result.v[1][2] ); 44 } 45 } 46 return 0; 47 } 48 // 求A^(n-1) 49 Matrix mtPow( int n ) 50 { 51 Matrix result,temp; 52 if(n == 1) return A ; 53 else 54 { 55 temp = mtPow( n/2 ) ; 56 result = mtMultiply( temp, temp ) ; 57 } 58 if(n % 2 != 0) 59 return mtMultiply(result , A) ; 60 else return result ; 61 } 62 63 Matrix mtMultiply( Matrix m, Matrix n ) 64 { 65 Matrix result ; 66 result.v[1][1] = (m.v[1][1]*n.v[1][1] + m.v[1][2]*n.v[2][1]) % 10000 ; 67 result.v[1][2] = (m.v[1][1]*n.v[1][2] + m.v[1][2]*n.v[2][2]) % 10000 ; 68 result.v[2][1] = (m.v[2][1]*n.v[1][1] + m.v[2][2]*n.v[2][1]) % 10000 ; 69 result.v[2][2] = (m.v[2][1]*n.v[1][2] + m.v[2][2]*n.v[2][2]) % 10000 ; 70 return result ; 71 }
事实上，如果A和Basic中如果有元素不会被使用到，我们应该尽可能地把它设成0，这样能加快计算的速度。
在排列组合中的应用：
题目链接：http://code.hdu.edu.cn/showproblem.php?pid=3519
HDU 3519 Lucky Coins Sequence
题目意思就不用我解释了吧。
这道题一看还以为是排列组合，但是一看数据量，10^9，无论是正常的DP还是排列组合，都是过不去的。
所以，总结了一下，大数DP可以考虑用矩阵快速幂。
可以先打表看一下规律：
显然不符合的数就是Fabbonaci的变形。
下面推一下DP递推式：
长度为 n 的 01 串一共有 2^n 种不同的排列方法 ,
设 f(n) 为长度是 n 的不包含连续 3 个或以上相同的 1 或 0 的 01 串 ,
则 f(1)=2,f(2)=4,f(3)=6,f(4)=10
当 n>4 的时候 , 分情况考虑 :
1、 如果是以 00 或者 11 结尾 , 则分别有 f(n-2)/2 种情况 , 加起来就是 f(n-2) 种 .
2、 如果是以 01 或者 10 结尾 , 则第 n 个字符要和第 n-1 个字符不一样 , 那么分别有 f(n-1)/2 种 , 加起来就是 f(n-1)
则统计起来就是 f(n)=f(n-1)+f(n-2), 题目要求的是包含连续三个相同的 0 或 1 串的串数 , 那就是用 a[n]=(2^n-f(n))%10007.
然而这样还不好求 , 先不看 %10007, 转换成递推公式是 a[n]=a[n-1]+a[n-2]+2^(n-2),
转换成矩阵 :
a[n] 1 1 1 a[n-1]
a[n-1] = 1 0 0 * a[n-2]
2^(n-1) 0 0 2 2^(n-2)
和上一题一样的方法，可以得到构造矩阵
| 1 1 1 | | 2 0 0 |
A = | 1 0 0 | Basic = | 0 0 0 | 则 F[n] = A^(n-3) * Basic ( n > 3 ) .
| 0 0 2 | | 4 0 0 |
n = 1,2,3 单独处理。
接下来的都不用多解释了吧，上代码：
1 #include<stdio.h> 2 #include<iostream> 3 using namespace std; 4 #include<queue> 5 #include<math.h> 6 #include<algorithm> 7 #include<string.h> 8 #include<stdlib.h> 9 #include<vector> 10 #include<set> 11 #include<map> 12 13 #define repA(p,q,i) for( int (i)=(p); (i)!=(q); ++(i) ) 14 #define repAE(p,q,i) for( int (i)=(p); (i)<=(q); ++(i) ) 15 #define repD(p,q,i) for( int (i)=(p); (i)!=(q); --(i) ) 16 #define repDE(p,q,i) for( int (i)=(p); (i)>=(q); --(i) ) 17 #define range 1010 18 19 20 struct Matrix{ 21 int v[4][4] ; 22 }; 23 24 Matrix basic,A; // 基础矩阵从F[3]开始 25 Matrix mtPow( int n ); 26 Matrix mtMultiply( Matrix m, Matrix n ); 27 28 int main() 29 { 30 // 基础矩阵从F[3]开始 31 memset( basic.v, 0, sizeof(basic.v) ) ; 32 basic.v[1][1]=2; 33 basic.v[2][1]=0; 34 basic.v[3][1]=4; 35 //递推矩阵A 36 memset( A.v, 0 , sizeof(basic.v) ) ; 37 A.v[1][1]=A.v[1][2]=A.v[1][3]=1; 38 A.v[2][1]=1; A.v[3][3]=2; 39 //solve 40 int n; 41 while( scanf("%d",&n) != EOF ) 42 { 43 if(n <=0 ) break; 44 else if(n==1 || n==2) printf("0\n"); 45 else if(n == 3) printf("2\n"); 46 else 47 { //注意矩阵相乘顺序 48 Matrix result = mtMultiply( mtPow(n-3) , basic ) ; 49 printf("%d\n", result.v[1][1] ); 50 } 51 } 52 return 0; 53 } 54 55 Matrix mtPow( int n ) 56 { 57 Matrix result,temp; 58 if(n == 1) return A ; 59 else 60 { 61 temp = mtPow(n/2) ; 62 result = mtMultiply( temp, temp ) ; 63 } 64 if(n % 2 != 0) 65 return mtMultiply(result , A) ; 66 else return result ; 67 } 68 69 Matrix mtMultiply( Matrix m, Matrix n ) 70 { 71 Matrix result ; 72 73 result.v[1][1] = (m.v[1][1]*n.v[1][1] + m.v[1][2]*n.v[2][1] + m.v[1][3]*n.v[3][1])%10007 ; 74 result.v[1][2] = (m.v[1][1]*n.v[1][2] + m.v[1][2]*n.v[2][2] + m.v[1][3]*n.v[3][2])%10007 ; 75 result.v[1][3] = (m.v[1][1]*n.v[1][3] + m.v[1][2]*n.v[2][3] + m.v[1][3]*n.v[3][3])%10007 ; 76 result.v[2][1] = (m.v[2][1]*n.v[1][1] + m.v[2][2]*n.v[2][1] + m.v[2][3]*n.v[3][1])%10007 ; 77 result.v[2][2] = (m.v[2][1]*n.v[1][2] + m.v[2][2]*n.v[2][2] + m.v[2][3]*n.v[3][2])%10007 ; 78 result.v[2][3] = (m.v[2][1]*n.v[1][3] + m.v[2][2]*n.v[2][3] + m.v[2][3]*n.v[3][3])%10007 ; 79 result.v[3][1] = (m.v[3][1]*n.v[1][1] + m.v[3][2]*n.v[2][1] + m.v[3][3]*n.v[3][1])%10007 ; 80 result.v[3][2] = (m.v[3][1]*n.v[1][2] + m.v[3][2]*n.v[2][2] + m.v[3][3]*n.v[3][2])%10007 ; 81 result.v[3][3] = (m.v[3][1]*n.v[1][3] + m.v[3][2]*n.v[2][3] + m.v[3][3]*n.v[3][3])%10007 ; 82 return result ; 83 }
所以大数DP，如果要用矩阵快速幂的话，总结起来就是打表，找规律，找递推式，构造矩阵。
关键在于DP递推式，其实找出递推式之后，基本上都是模板了。
而HDU 3519 更为一般的解法是：
