• POJ 3233 Matrix Power Series (矩阵 && 求和 && 线性变换)


    题目大意:给定n阶方阵A,计算S = A^1+A^2+……+A^k (mod m)  (k <= 10^9)   我们可以根据矩阵快速幂在O(n^3log(k))的时间里算出A^k,但我们也不能一个一个算再加起来啊,那样铁定超时……   ①二分 这种方法我觉得与秦九韶算法计算多项式的思路类似,都是找出重复因子而减少多项式计算的次数.当然我们这个多项式比较特殊,所以有比秦九韶算法更好的思路: 当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   这种方法可能出现递归过深的问题,所以要设置好占空间小心栈溢出,而且速度也偏慢
    11217405 AbandonZHANG 3233 Accepted 3532K 704MS G++ 2444B 2013-01-28 20:53:36
     
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #define MID(x,y) ((x+y)>>1)
    #pragma comment(linker, "/STACK:102400000,102400000")
    using namespace std;
    typedef long long LL;
    
    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)   这样要比①好多了,又省空间又省时间。
    11217490 AbandonZHANG 3233 Accepted 1016K 282MS G++ 2962B 2013-01-28 21:25:37
     
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #include 
    #define MID(x,y) ((x+y)>>1)
    #pragma comment(linker, "/STACK:102400000,102400000")
    using namespace std;
    typedef long long LL;
    
    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;
    }
    
     
    举杯独醉,饮罢飞雪,茫然又一年岁。 ------AbandonZHANG
  • 相关阅读:
    yum install mysql.i686
    firefox无法浏览flash的解决方案
    vb.net如何打开指定文件
    XML文件操作的简单类
    window server 安装与卸载
    常用的sql语句
    with进行递归表
    常用js
    调用Google的自动翻译
    MySQl 总结知识
  • 原文地址:https://www.cnblogs.com/AbandonZHANG/p/4113992.html
Copyright © 2020-2023  润新知