POJ 3233 -- Matrix Power Series
Description
Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.
Input
The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.
Output
Output the elements of S modulo m in the same way as A is given.
Sample Input
2 2 4 0 1 1 1
Sample Output
1 2 2 3
http://www.abandonzhang.com/archives/484
题目给定一个n*n的矩阵A，要求(A+A^2+....A^K)%m后的矩阵
①二分
这种方法我觉得与秦九韶算法计算多项式的思路类似，都是找出重复因子而减少多项式计算的次数.当然我们这个多项式比较特殊，所以有比秦九韶算法更好的思路：
当k为偶数时：
S(k) = A^1+A^2+A^3 + A^(k/2) + A^(k/2+1) + …… + A^(k/2+k/2)
= (A^(k/2) + E) (A^1 + A^2 + A^3 …… + A^(k/2) )
=(A^(k/2) + E) * S(k/2)
当k为奇数时：
S(k) = S(k-1) * A^k
const int MAX = 30;
struct Mat{
int row, col;
LL mat[MAX][MAX];
};
//initialize square matrix to unit matrix
Mat unit(int n){
Mat A;
A.row = A.col = n;
memset(A.mat, 0, sizeof(A.mat));
for (int i = 0; i < n; i ++)
A.mat[i][i] = 1;
return A;
}
//return T(A)
Mat transpose(Mat A){
Mat T;
T.row = A.col;
T.col = A.row;
for (int i = 0; i < A.col; i ++)
for (int j = 0; j < A.row; j ++)
T.mat[i][j] = A.mat[j][i];
return T;
}
//return (A+B)%mod
Mat add(Mat A, Mat B, int mod){
Mat C = A;
for (int i = 0; i < C.row; i ++)
for (int j = 0; j < C.col; j ++)
C.mat[i][j] = (A.mat[i][j] + B.mat[i][j]) % mod;
return C;
}
//return A*B%mod
Mat mul(Mat A, Mat B, int mod){
Mat C;
C.row = A.row;
C.col = B.col;
for (int i = 0; i < A.row; i ++){
for (int j = 0; j < B.col; j ++){
C.mat[i][j] = 0;
for (int k = 0; k < A.col; k ++)
//注意这里要保证乘法不溢出，否则还需要设计特殊的乘法模
C.mat[i][j] += A.mat[i][k] * B.mat[k][j];
C.mat[i][j] %= mod;
}
}
return C;
};
//return A^n%mod
Mat exp_mod(Mat A, int n, int mod){
Mat res = unit(A.row);
while(n){
if (n & 1){
res = mul(res, A, mod);
}
A = mul(A, A, mod);
n >>= 1;
}
return res;
}
Mat A;
int n, k, mod;
Mat S(int k){
if (k == 1){
return A;
}
if (k % 2 == 0){
return mul(add(exp_mod(A, k/2, mod), unit(n), mod), S(k/2), mod);
}
else{
return add(S(k-1), exp_mod(A, k, mod), mod);
}
}
int main(){
cin >> n >> k >> mod;
A.row = A.col = n;
for (int i = 0; i < n; i ++)
for (int j = 0; j < n; j ++)
cin >> A.mat[i][j];
Mat ans = S(k);
for (int i = 0; i < ans.row; i ++){
for (int j = 0; j < ans.col - 1; j ++)
cout << ans.mat[i][j] << " ";
cout << ans.mat[i][ans.col-1] << endl;
}
return 0;
}
②线性变换(很赞的一种方法呐~！涨姿势~)
我们考虑递推，即从s(k-1)到s(k)的线性变换。
首先一维线性变换显然是出不来的，就像A^1+A^2+A^3+……+A^k = d * (A^1+A^2+A^3+……+A^(k-1) ) 明显不行……
那么我们再加一维就显然了：
A^1+A^2+A^3+……+A^k =( A^1+A^2+A^3+……+A^(k-1) ) + A^k
A^(k+1) = 0 * ( A^1+A^2+A^3+……+A^(k-1) ) + A * A^k. //这一行的变换是作为辅助，补全二维的线性变换.
那么线性变换矩阵B显然就是：
1 1
0 1
一般化就有
s(k) 1 1 s(k-1)
A^k+1 = 0 1 * A^k
即P(k) = B^(n-1) * P(1)
const int MAX = 60;
struct Mat{
int row, col;
LL mat[MAX][MAX];
};
//initialize square matrix to unit matrix
Mat unit(int n){
Mat A;
A.row = A.col = n;
memset(A.mat, 0, sizeof(A.mat));
for (int i = 0; i < n; i ++)
A.mat[i][i] = 1;
return A;
}
//return T(A)
Mat transpose(Mat A){
Mat T;
T.row = A.col;
T.col = A.row;
for (int i = 0; i < A.col; i ++)
for (int j = 0; j < A.row; j ++)
T.mat[i][j] = A.mat[j][i];
return T;
}
//return (A+B)%mod
Mat add(Mat A, Mat B, int mod){
Mat C = A;
for (int i = 0; i < C.row; i ++)
for (int j = 0; j < C.col; j ++)
C.mat[i][j] = (A.mat[i][j] + B.mat[i][j]) % mod;
return C;
}
//return A*B%mod
Mat mul(Mat A, Mat B, int mod){
Mat C;
C.row = A.row;
C.col = B.col;
for (int i = 0; i < A.row; i ++){
for (int j = 0; j < B.col; j ++){
C.mat[i][j] = 0;
for (int k = 0; k < A.col; k ++)
//注意这里要保证乘法不溢出，否则还需要设计特殊的乘法模
C.mat[i][j] += A.mat[i][k] * B.mat[k][j];
C.mat[i][j] %= mod;
}
}
return C;
};
//return A^n%mod
Mat exp_mod(Mat A, int n, int mod){
Mat res = unit(A.row);
while(n){
if (n & 1){
res = mul(res, A, mod);
}
A = mul(A, A, mod);
n >>= 1;
}
return res;
}
Mat A;
int n, k, mod;
int main(){
cin >> n >> k >> mod;
A.row = A.col = n;
for (int i = 0; i < n; i ++)
for (int j = 0; j < n; j ++)
cin >> A.mat[i][j];
Mat res;
res.row = 2 * n;
res.col = n;
for (int i = 0; i < n; i ++)
for (int j = 0; j < n; j ++)
res.mat[i][j] = A.mat[i][j];
for (int i = n; i < 2 * n; i ++)
for (int j = 0; j < n; j ++)
res.mat[i][j] = A.mat[i-n][j];
Mat E = unit(n);
Mat B;
memset(B.mat, 0, sizeof(B.mat));
B.row = B.col = 2 * n;
for (int i = 0; i < n; i ++){
for (int j = 0; j < n; j ++)
B.mat[i][j] = E.mat[i][j];
for (int j = n; j < 2 * n; j ++)
B.mat[i][j] = A.mat[i][j-n];
}
for (int i = n; i < 2 * n; i ++){
for (int j = n; j < 2 * n; j ++)
B.mat[i][j] = A.mat[i-n][j-n];
}
B = exp_mod(B, k-1, mod);
res = mul(B, res, mod);
for (int i = 0; i < n; i ++){
for (int j = 0; j < n - 1; j ++)
cout << res.mat[i][j] << " ";
cout << res.mat[i][n-1] << endl;
}
return 0;
}
对于任意的A^x，我们都能够利用矩阵快速幂求出，但是我们现在要求的是和。
仔细观察整个式子，那么我们可以对原式进行变形
如果k为偶数，那么(A+A^2+....A^K) = (A+...+A^K/2)+A^K/2*(A+...+A^K/2)
如果k为奇数，那么(A+A^2+....A^K) = (A+...+A^K/2)+A^K/2*(A+...+A^K/2)+A^k
3 那么对于上面的式子的变形，就是二分的思想，那么我们可以利用二分来求和，然后对于单个的矩阵的x次方我们利用快速幂
思路：矩阵快速幂。首先我们知道 A^x 可以用矩阵快速幂求出来（具体可见poj 3070）。其次可以对k进行二分，每次将规模减半，分k为奇偶两种情况，如当k = 6和k = 7时有：
ps：对矩阵定义成结构体Matrix，求S时用递归，程序会比较直观，好写一点。当然定义成数组，然后再进行一些预处理，效率会更高些。
- const int N = 30;
- int n , k , MOD;
- struct Matrix{
- int mat[N][N];
- Matrix operator*(const Matrix& m)const{
- Matrix tmp;
- for(int i = 0 ; i < n ; i++){
- for(int j = 0 ; j < n ; j++){
- tmp.mat[i][j] = 0;
- for(int k = 0 ; k < n ; k++)
- tmp.mat[i][j] += mat[i][k]*m.mat[k][j]%MOD;
- tmp.mat[i][j] %= MOD;
- }
- }
- return tmp;
- }
- Matrix operator+(const Matrix& m)const{
- Matrix tmp;
- for(int i = 0 ; i < n ; i++)
- for(int j = 0 ; j < n ; j++)
- tmp.mat[i][j] = (mat[i][j]+m.mat[i][j])%MOD;
- return tmp;
- }
- };
- Matrix Pow(Matrix m , int t){
- Matrix ans;
- memset(ans.mat , 0 , sizeof(ans.mat));
- for(int i = 0 ; i < n ; i++)
- ans.mat[i][i] = 1;
- while(t){
- if(t&1)
- ans = ans*m;
- t >>= 1;
- m = m*m;
- }
- return ans;
- }
- Matrix solve(Matrix m , int t){
- Matrix A;
- memset(A.mat , 0 , sizeof(A.mat));
- for(int i = 0 ; i < n ; i++)
- A.mat[i][i] = 1;
- if(t == 1)
- return m;
- if(t&1)
- return (Pow(m,t>>1)+A)*solve(m,t>>1)+Pow(m , t);
- else
- return (Pow(m,t>>1)+A)*solve(m,t>>1);
- }
- void output(Matrix ans){
- for(int i = 0 ; i < n ; i++){
- printf("%d" , ans.mat[i][0]%MOD);
- for(int j = 1 ; j < n ; j++)
- printf(" %d" , ans.mat[i][j]%MOD);
- puts("");
- }
- }
- int main(){
- Matrix m , ans;
- while(scanf("%d%d%d" , &n , &k , &MOD) != EOF){
- for(int i = 0 ; i < n ; i++)
- for(int j = 0 ; j < n ; j++)
- scanf("%d" , &m.mat[i][j]);
- ans = solve(m , k);
- output(ans);
- }
- return 0;
- }
先把形式写成 Sk = A + A2 + A3 + … + Ak, 然后容易看出
Sk = Sk/2 × Ak/2 + Sk/2+ (Ak if odd k)
这样的解法用矩阵快速幂已经可以通过，但是每次还要矩阵加法。
因为 Ak = A × Ak-1，如果把 A 看成常数，然后 Sk 也用 Sk-1 和 Ak-1 表示出来，就可以把数列的递推式用矩阵表示。
于是 Sk= A × Ak-1 + Sk-1，所以
这样只要算一个快速幂就好了，矩阵规模被扩大了两倍，但是根据这个矩阵的特点算乘法的时候可以做些优化
