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
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/AbandonZHANG/p/4114210.html