POJ 3233 -- Matrix Power Series

POJ 3233 -- Matrix Power Series

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 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;
}

http://blog.csdn.net/chenguolinblog/article/details/10454943

对于任意的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时有：

      k = 6 有： S(6) = (1 + A^3) * (A + A^2 + A^3) = (1 + A^3) * S(3)。

      k = 7 有： S(7) = A + (A + A^4) * (A + A^2 + A^3) = A + (A + A^4) * S(3)。

ps：对矩阵定义成结构体Matrix，求S时用递归，程序会比较直观，好写一点。当然定义成数组，然后再进行一些预处理，效率会更高些。

1. const int N = 30;  
2.   
3. int n , k , MOD;  
4. struct Matrix{  
5.     int mat[N][N];  
6.     Matrix operator*(const Matrix& m)const{  
7.         Matrix tmp;  
8.         for(int i = 0 ; i < n ; i++){  
9.             for(int j = 0 ; j < n ; j++){  
10.                 tmp.mat[i][j] = 0;  
11.                 for(int k = 0 ; k < n ; k++)  
12.                     tmp.mat[i][j] += mat[i][k]*m.mat[k][j]%MOD;  
13.                 tmp.mat[i][j] %= MOD;  
14.             }  
15.         }  
16.         return tmp;  
17.     }  
18.     Matrix operator+(const Matrix& m)const{  
19.         Matrix tmp;  
20.         for(int i = 0 ; i < n ; i++)  
21.             for(int j = 0 ; j < n ; j++)  
22.                 tmp.mat[i][j] = (mat[i][j]+m.mat[i][j])%MOD;  
23.         return tmp;  
24.     }  
25. };  
26.   
27. Matrix Pow(Matrix m , int t){  
28.     Matrix ans;  
29.     memset(ans.mat , 0 , sizeof(ans.mat));  
30.     for(int i = 0 ; i < n ; i++)  
31.         ans.mat[i][i] = 1;  
32.     while(t){  
33.         if(t&1)  
34.             ans = ans*m;  
35.         t >>= 1;  
36.         m = m*m;  
37.     }  
38.     return ans;  
39. }  
40.   
41. Matrix solve(Matrix m , int t){  
42.     Matrix A;  
43.     memset(A.mat , 0 , sizeof(A.mat));  
44.     for(int i = 0 ; i < n ; i++)  
45.         A.mat[i][i] = 1;  
46.     if(t == 1)  
47.         return m;  
48.     if(t&1)  
49.         return (Pow(m,t>>1)+A)*solve(m,t>>1)+Pow(m , t);         
50.     else  
51.         return (Pow(m,t>>1)+A)*solve(m,t>>1);         
52. }  
53.   
54. void output(Matrix ans){  
55.     for(int i = 0 ; i < n ; i++){  
56.         printf("%d" , ans.mat[i][0]%MOD);  
57.         for(int j = 1 ; j < n ; j++)  
58.             printf(" %d" , ans.mat[i][j]%MOD);  
59.         puts("");  
60.     }  
61. }  
62.   
63. int main(){  
64.     Matrix m , ans;  
65.     while(scanf("%d%d%d" , &n , &k , &MOD) != EOF){  
66.          for(int i = 0 ; i < n ; i++)  
67.             for(int j = 0 ; j < n ; j++)  
68.                 scanf("%d" , &m.mat[i][j]);  
69.          ans = solve(m , k);  
70.          output(ans);  
71.     }  
72.     return 0;  
73. }  

http://xufz.wordpress.com/2013/12/21/poj-3233-matrix-power-series/

先把形式写成 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，所以

这样只要算一个快速幂就好了，矩阵规模被扩大了两倍，但是根据这个矩阵的特点算乘法的时候可以做些优化

Read full article from 3233 -- Matrix Power Series