算法习题——矩阵
Q1(poj 3070):
求斐波那契数列的第n个，n最大取到1000000000,。
分析：求这种较大递推数列的一般方法使用矩阵快速幂的，这里题目直接给出了矩阵形式，就不需要进行友矩阵(A)的构造了，也不需要进行最后一次矩阵和向量的相乘，直接初始化矩阵规模，进行快速幂即可。
不过这个结论倒是可以记住，以后就有了logn求斐波那契数列的方法：
参考代码如下：
#include<iostream> #include<string> #include<cstring> using namespace std; const int maxn = 20; typedef long long Matrix[maxn][maxn]; typedef long long Vector[maxn]; int sz , mod; void matrix_mul(Matrix A , Matrix B , Matrix res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_pow(Matrix A , int n , Matrix res)//矩阵快速幂O(logn) { Matrix a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } void Transform(Vector d , Matrix A , Vector res)//计算完A^(n-d)之后，这个函数计算A*d = res { Vector r; memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz;j++) r[i] = (r[i] + d[j] * A[i][j])%mod; memcpy(res , r , sizeof(r)); } int main() { int n; while(cin >>n) { if(n == -1) break; if(n == 0) {cout<<0<<endl;continue;} mod = 10000; sz = 2; Vector a , f; Matrix A; memset(A , 0 , sizeof(A)); A[0][0] = 1; A[0][1] = 1; A[1][0] = 1; A[1][1] = 0; matrix_pow(A , n - 1 , A); cout << A[0][0] << endl; } return 0; }
Q2(hdu 5950):
一个数列的递推关系式是f[i] = f[i-1] + 2f[i-2] + i^4，现在给出n(最大可取1000000000)，计算f[n]%mod的结果。
分析：这道题的核心方法就是上面我们提到的时间复杂度为O(logn)求递推关系的矩阵快速幂。而这个方法的关键就是快速幂的底数A，需要一定的推导技巧。
#include<iostream> #include<string> #include<cstring> using namespace std; const int maxn = 20; typedef long long Matrix[maxn][maxn]; typedef long long Vector[maxn]; int sz ; long long mod = 2147493647; void matrix_mul(Matrix A , Matrix B , Matrix res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_pow(Matrix A , int n , Matrix res)//矩阵快速幂O(logn) { Matrix a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } void Transform(Vector d , Matrix A , Vector res)//计算完A^(n-d)之后，这个函数计算A*d = res { Vector r; memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz;j++) r[i] = (r[i] + d[j] * A[i][j])%mod; memcpy(res , r , sizeof(r)); } int main() { long long n , a , b; int t; cin>>t; while(t--) { cin>>n>>a>>b; sz = 7; //矩阵A的规模是sz Vector f; if(n == 1) {cout<<a<<endl;continue;} if(n == 2) {cout<<b<<endl;continue;} f[0] = b;f[1] = a; f[2] = 16;f[3] = 8;f[4] = 4;f[5] = 2;f[6] = 1; Matrix A; memset(A , 0 , sizeof(A));//注意每次矩阵A的参数都要初始化. A[0][0] = 1;A[0][1] = 2;A[0][2] = 1;A[0][3] = 4;A[0][4] = 6;A[0][5] = 4;A[0][6] = 1; A[1][0] = 1; A[2][2] = 1;A[2][3] = 4;A[2][4] = 6;A[2][5] = 4;A[2][6] = 1; A[3][3] = 1;A[3][4] = 3;A[3][5] = 3;A[3][6] = 1; A[4][4] = 1;A[4][5] = 2;A[4][6] = 1; A[5][5] = 1;A[5][6] = 1; A[6][6] = 1; matrix_pow(A , n - 2 , A); Transform(f , A , f); cout << f[0] << endl; } return 0; }
Q3(uva 10689):
这道题目是矩阵快速幂求斐波那契数列的推广形式，给出f[0],f[1]的初值，递推关系式斐波那契数列的递推关系f[n]=f[n-1]+f[n-2].
分析：由于递推关系和斐波那契数列的递推关系一致，友矩阵Q是一样的.
#include<iostream> #include<string> #include<cstring> #include<cmath> using namespace std; const int maxn = 20; typedef long long Matrix[maxn][maxn]; typedef long long Vector[maxn]; int sz , mod; void matrix_mul(Matrix A , Matrix B , Matrix res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_pow(Matrix A , int n , Matrix res)//矩阵快速幂O(logn) { Matrix a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } void Transform(Vector d , Matrix A , Vector res)//计算完A^(n-d)之后，这个函数计算A*d = res { Vector r; memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz;j++) r[i] = (r[i] + d[j] * A[i][j])%mod; memcpy(res , r , sizeof(r)); } int main() { int a , b , m , n; int t; cin>>t; while(t--) { cin >> a >> b >> n >>m; mod = 10; sz = 2; for(int i = 1;i < m;i++) mod *= 10; Vector f; Matrix A; f[1] = a; f[0] = b; //列向量[f(d) , f(d-1) , f(d-2) , ... , f(1)] memset(A , 0 , sizeof(A));//得到配出来的常数矩阵A A[0][0] = 1;A[0][1] = 1; A[1][0] = 1;A[1][1] = 0; matrix_pow(A , n - 1 , A); Transform(f , A , f); cout << f[0] << endl; } return 0; }
Q4(hdu4565):
Q5(hdu2256)：
#include<iostream> #include<string> #include<cstring> #include<cstdio> using namespace std; const int maxn = 20; typedef long long Matrix[maxn][maxn]; typedef long long Vector[maxn]; int sz , mod; void matrix_mul(Matrix A , Matrix B , Matrix res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_pow(Matrix A , int n , Matrix res)//矩阵快速幂O(logn) { Matrix a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } void Transform(Vector d , Matrix A , Vector res)//计算完A^(n-d)之后，这个函数计算A*d = res { Vector r; memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz;j++) r[i] = (r[i] + d[j] * A[i][j])%mod; memcpy(res , r , sizeof(r)); } int main(){ int n , T; scanf("%d" , &T); while(T--){ scanf("%d" , &n); mod = 1024;//整个数值在m的剩余系下 sz = 2;//矩阵A的规模是sz if(n == 0) {printf("1\n");continue;} Vector a , f; Matrix A; //A是矩阵,a是递推公式的系数向量 f是数列的初值向量 f[0] = 1; f[1] = 0; //列向量[f(d) , f(d-1) , f(d-2) , ... , f(1)] //数学概念上向量从1开始，但是存储上f[]数组是从0开始的，这点一定要注意!!逻辑高位存储在存储结构的地位，这里非常容易错！ memset(A , 0 , sizeof(A));//得到配出来的常数矩阵A A[0][0] = 5;A[0][1] = 12; A[1][0] = 2;A[1][1] = 5; matrix_pow(A , n , A); Transform(f , A , f); printf("%d\n" ,(2 * f[0] - 1) %mod); } return 0; }
Q6(hdu5451):
#include<iostream> #include<string> #include<cstring> #include<cstdio> using namespace std; const int maxn = 20; typedef long long Matrix[maxn][maxn]; typedef long long Vector[maxn]; typedef long long LL; int sz , mod; void matrix_mul(Matrix A , Matrix B , Matrix res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_pow(Matrix A , int n , Matrix res)//矩阵快速幂O(logn) { Matrix a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } void Transform(Vector d , Matrix A , Vector res)//计算完A^(n-d)之后，这个函数计算A*d = res { Vector r; memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz;j++) r[i] = (r[i] + d[j] * A[i][j])%mod; memcpy(res , r , sizeof(r)); } LL Quick_Mod(LL a, LL b, int m){ LL res = 1,term = a % m; while(b) { if(b & 1) res = (res * term) % m; term = (term * term) % m; b >>= 1; } return res%m; } int main(){ int T; LL x; scanf("%d" , &T); for(int kase = 1;kase <= T;kase++){ scanf("%lld %d" , &x , &mod); sz = 2;//矩阵A的规模是sz Vector a , f; Matrix A; //A是矩阵,a是递推公式的系数向量 f是数列的初值向量 f[0] = 1; f[1] = 0; x = Quick_Mod(2 , x , mod * mod - 1) + 1; if(x == 0) {printf("1\n");continue;} //列向量[f(d) , f(d-1) , f(d-2) , ... , f(1)] //数学概念上向量从1开始，但是存储上f[]数组是从0开始的，这点一定要注意!!逻辑高位存储在存储结构的地位，这里非常容易错！ memset(A , 0 , sizeof(A));//得到配出来的常数矩阵A A[0][0] = 5;A[0][1] = 12; A[1][0] = 2;A[1][1] = 5; matrix_pow(A , x , A); Transform(f , A , f); printf("Case #%d: %d\n" ,kase , (2 * f[0] - 1) %mod); } return 0; }
Q7(hdu6198):
基于斐波那契数列，给出整数k，让你计算k个允许重复的斐波那契数列，不能组成的最小整数是多少。
分析：通过手动算找规律，编程打表找规律也好，找到规律问题的本质是让我们求fib(2n + 3)-1.然后简单一发矩阵快速幂即可。
Q8(hdu6185):
用1x2的格子覆盖4xn的格子，问有多少种不同的方案，结果对1e9+7取模。
分析：这道题目算是很综合了。纯粹基于覆盖的机理找递推关系难度较大，可以通过dfs找到小数据的表，我们根据n的取值，也能够猜测这道题目应该是有线性递推的关系的。这里在拿到了前十项的数据之后，找规律有个小套路，最基本的是靠观察，但是对于含有2项以上的线性递推关系，这个方法就不是很好了。这里我们依次假设题目有如下的线性递推关系：
(1)f(n) = a*f(n-1)+b*f(n-2)
(2)f(n) = a*f(n-1)+b*f(n-2)+c*f(n-3)
(3)f(n)=a*f(n-1)+b*f(n-2)+c*f(n-3)+d*f(n-4)
...
然后基于原始数据，构造线性方程组再用高斯消元求解a、b等参数进行验证，当得到系数均为整数的时候，即题目满足的线性递推关系。
找到了线性递推关系之后，用矩阵快速幂跑一发即可。
Q9(hdu4965):
#include<iostream> #include<string> #include<cstring> #include<cstdio> using namespace std; const int maxn = 1000 + 1; typedef long long Matrix1[maxn][6]; typedef long long Matrix2[6][maxn]; typedef long long Matrix3[6][6]; int sz , mod; void matrix_mul(Matrix3 A , Matrix3 B , Matrix3 res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix3 C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_mul_p1(Matrix2 B , Matrix1 A, int sz1 , int sz2, int sz3 , Matrix3 res){//sz1*sz2的矩阵与sz2*sz3的矩阵做乘法 Matrix3 C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz1;i++) for(int j = 0;j < sz3;j++) for(int k = 0;k < sz2;k++) C[i][j] = (C[i][j] + B[i][k]*A[k][j]); memcpy(res , C , sizeof(C)); } void matrix_mul_p2(Matrix1 A , Matrix3 C, int sz1 , int sz2, int sz3 , Matrix1 res){//sz1*sz2的矩阵与sz2*sz3的矩阵做乘法 Matrix1 temp; memset(temp , 0 , sizeof(temp)); for(int i = 0;i < sz1;i++) for(int j = 0;j < sz3;j++) for(int k = 0;k < sz2;k++) temp[i][j] = (temp[i][j] + A[i][k]*C[k][j])%mod; memcpy(res , temp , sizeof(temp)); } int matrix_mul_p3(Matrix1 A , Matrix2 B, int sz1 , int sz2, int sz3){//sz1*sz2的矩阵与sz2*sz3的矩阵做乘法 int temp ,ans; temp = 0;ans = 0; for(int i = 0;i < sz1;i++){ for(int j = 0;j < sz3;j++){ temp = 0; for(int k = 0;k < sz2;k++) temp = (temp + A[i][k]*B[k][j])%mod; ans += temp; } } return ans; } void matrix_pow(Matrix3 A , int n , Matrix3 res)//矩阵快速幂O(logn) { Matrix3 a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } int main(){ int n , k; while(scanf("%d%d" , &n , &k)){ Matrix1 A; Matrix2 B; Matrix3 C; mod = 6; for(int i = 0;i < n;i++) for(int j = 0;j < k;j++) scanf("%lld" , &A[i][j]); for(int i = 0;i < k;i++) for(int j = 0;j < n;j++) scanf("%lld" , &B[i][j]); matrix_mul_p1(B , A , k , n , k , C); sz = k; matrix_pow(C , n*n - 1, C); //下面要做n*k的A * k*k的C * k*n的B matrix_mul_p2(A , C , n , k , k , A); int ans = matrix_mul_p3(A , B , n , k , n); printf("%d\n" , ans); } return 0; } //lrj的快速幂模板迷之超时，下面贴一个ac代码 /* #include <cstdio> #include <cstring> #include <cmath> #include <cstdlib> #include <iostream> #include <algorithm> #include <queue> #include <map> #include <set> #include <vector> #include <cctype> using namespace std; const int mod=6; struct matrix{ int f[6][6]; }; int A[1001][6],B[6][1001],C[1001][6],D[1001][1001]; matrix mul(matrix a,matrix b,int n) { matrix c; memset(c.f,0,sizeof(c.f)); int i,j,k; for(i=0;i<n;i++) { for(j=0;j<n;j++) { for(k=0;k<n;k++) { c.f[i][j]+=a.f[i][k]*b.f[k][j]; } c.f[i][j]%=mod; } } return c; } matrix pow_mod(matrix a,int b,int n) { matrix s; memset(s.f,0,sizeof(s.f)); for(int i=0;i<n;i++)s.f[i][i]=1; while(b) { if(b&1)s=mul(s,a,n); a=mul(a,a,n); b=b>>1; } return s; } int main() { int n,K; while(scanf("%d%d",&n,&K)!=EOF) { if(n==0&&K==0)break; int i,j,k; for(i=0;i<n;i++) for(j=0;j<K;j++) scanf("%d",&A[i][j]); for(i=0;i<K;i++) for(j=0;j<n;j++) scanf("%d",&B[i][j]); matrix e,g; memset(e.f,0,sizeof(e.f)); for(i=0;i<K;i++) { for(j=0;j<K;j++) { for(k=0;k<n;k++) e.f[i][j]+=B[i][k]*A[k][j]; e.f[i][j]%mod; } } e=pow_mod(e,n*n-1,K); for(int i = 0;i < K;i++){ for(int j =0;j < K;j++) printf("%d " , e.f[i][j]); printf("\n"); } memset(C,0,sizeof(C)); for(i=0;i<n;i++) { for(j=0;j<K;j++) { for(k=0;k<K;k++) C[i][j]+=A[i][k]*e.f[k][j]; C[i][j]%=mod; } } for(int i = 0;i < n;i++){ for(int j =0;j < K;j++) printf("%d " , C[i][j]); printf("\n"); } int ans=0; memset(D,0,sizeof(D)); for(i=0;i<n;i++) { for(j=0;j<n;j++) { for(k=0;k<K;k++) D[i][j]+=C[i][k]*B[k][j]; ans+=D[i][j]%mod; } } printf("%d\n",ans); } return 0; } */
Q10(hdu4990):
利用题目给出的代码所呈现的机制，输入n、m输出对应的答案。
分析：这是一个典型的优化代码题目。能够看到m是mod运算使用的，可以不考虑。然后对于不同的n，我们通过题目给出的程序可以得到n取[0,11]的答案。
a[n]:0 1 2 5 10 21 42 85 170 341 682 1365.
题目其实就是快速求解这个数列.根据观察得，这个数列满足f[n] = f[n-1]+2f[n-2]+1的递推关系。然后矩阵快速幂构造一发即可。
#include<cstdio> #include<string> #include<cstring> using namespace std; const int maxn = 20; typedef long long Matrix[maxn][maxn]; typedef long long Vector[maxn]; int sz , mod; void matrix_mul(Matrix A , Matrix B , Matrix res)//矩阵A*B = C , O(n^3).矩阵均为同阶方阵 { Matrix C; memset(C , 0 , sizeof(C)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz; j++) for(int k = 0;k < sz; k++) C[i][j] = (C[i][j] + A[i][k] * B[k][j])%mod; memcpy(res , C , sizeof(C)); } void matrix_pow(Matrix A , int n , Matrix res)//矩阵快速幂O(logn) { Matrix a , r; memcpy(a , A , sizeof(a)); memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) r[i][i] = 1;//单位矩阵 while(n){ if(n&1) matrix_mul(r , a , r); n >>= 1; matrix_mul(a , a , a); } memcpy(res , r , sizeof(r)); } void Transform(Vector d , Matrix A , Vector res)//计算完A^(n-d)之后，这个函数计算A*d = res { Vector r; memset(r , 0 , sizeof(r)); for(int i = 0;i < sz;i++) for(int j = 0;j < sz;j++) r[i] = (r[i] + d[j] * A[i][j])%mod; memcpy(res , r , sizeof(r)); } int main(){ int n , m; while(~scanf("%d%d" , &n ,&m)){ mod = m;sz = 4; Vector a , f; Matrix A; f[0] = 2; f[1] = 1; f[2] = 0; f[3] = 1; memset(A , 0 , sizeof(A)); A[0][0] = 1;A[0][1] = 2;A[0][2] = 0;A[0][3] = 1; A[1][0] = 1;A[1][1] = 0;A[1][2] = 0;A[1][3] = 0; A[2][0] = 0;A[2][1] = 1;A[2][2] = 0;A[2][3] = 0; A[3][0] = 0;A[3][1] = 0;A[3][2] = 0;A[3][3] = 1; if(n ==0) {printf("0\n"); continue;} else if(n == 1) {printf("%d\n" , 1 % mod); continue;} else if(n == 2) {printf("%d\n" , 2 % mod); continue;} matrix_pow(A , n - 2 , A); Transform(f , A , f); printf("%d\n" , f[0]); } return 0; }
当然还有比较傻瓜的做法是，根据手头的数据，推测线性递推关系然后构造线性方程组再用高斯消元(上文一个题目用到过)。这里值得注意的是，推测线性递推关系中，应该考虑常数项的，这样相对上面(Q8:hdu6185)，实际就是多了一个自变量。
高斯消元的代码如下：
#include <iostream> #include <cstdio> #include <cmath> #include <cstring> #include <algorithm> #include <cstring> #include <vector> using namespace std; const int maxn=105; typedef double Matrix[maxn][maxn]; Matrix A,S; //n是方程的个数 void gauss(Matrix A,int n){ int i,j,k,r; for(int i=0; i<n; i++) { r=i; for( j=i+1; j<n; j++) if(fabs(A[j][i])>fabs(A[r][i]))r=j; if(r!=i) for(j=0; j<=n; j++)swap(A[r][j],A[i][j]); for(k=i+1; k<n; k++) { double f=A[k][i]/A[i][i]; for(j=i; j<=n; j++) A[k][j]-=f*A[i][j]; } } for(i=n-1; i>=0; i--) { for(j=i+1; j<n; j++) A[i][n]-=A[j][n]*A[i][j]; A[i][n]/=A[i][i]; } } int main() { memset(A,0,sizeof(A)); A[0][0]=1,A[0][1]=0,A[0][2]=1; A[1][0]=2,A[1][1]=1,A[1][2]=1; A[2][0]=5,A[2][1]=2,A[2][2]=1; //对于线性方程组Ax = b 代码中的A是[A,b] A[0][3]=2,A[1][3]=5,A[2][3]=10; gauss(A,3); //经过该代码实现的高斯消元，解向量x是A矩阵的最后一列，即替代了原来的系数向量 for(int i=0; i<3; i++) printf("%8.2f\n",A[i][3]); return 0; }
可求方程组的解集是[1,2,1]，对应着我们构造线性方程组的方式，可得到递推关系：f[n] = f[n-1] + 2f[n-2]+1.
