UVa 11149 Power of Matrix 矩阵快速幂-优快云博客

题意：

给出一个\(n \times n\)的矩阵\(A\)，求\(A+A^2+A^3+ \cdots + A^k\)。

分析：

这题是有\(k=0\)的情况，我们一开始先特判一下，直接输出单位矩阵\(E\)。
下面讨论\(k > 0\)的情况：

方法一

设答案为\(S_k(k > 0)\)
把矩阵增广一下
\(\begin{bmatrix} A & O \\ E & E \end{bmatrix} \begin{bmatrix} A^n\\ S_{n-1} \end{bmatrix} = \begin{bmatrix} A^{n+1}\\ S_n \end{bmatrix}\)
\(E\)表示单位矩阵，\(O\)是全为零的矩阵。
\({\begin{bmatrix} A & O \\ E & E \end{bmatrix}}^k \begin{bmatrix} A\\ O \end{bmatrix} = \begin{bmatrix} A^{k+1}\\ S_k \end{bmatrix}\)

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int maxn = 40;
const int MOD = 10;

int mul(int a, int b) { return a * b % MOD; }

void add(int& a, int b) { a += b; if(a >= MOD) a -= MOD; }

int n, k, sz;
int a[maxn][maxn], b[maxn][maxn];

struct Matrix
{
    int a[maxn * 2][maxn * 2];

    Matrix() { memset(a, 0, sizeof(a)); }

    Matrix operator * (const Matrix& t) const {
        Matrix ans;
        for(int i = 0; i < sz; i++)
            for(int j = 0; j < sz; j++)
                for(int k = 0; k < sz; k++)
                    add(ans.a[i][j], mul(a[i][k], t.a[k][j]));
        return ans;
    }

    void output() {
        for(int i = 0; i < sz; i++) {
            for(int j = 0; j < sz - 1; j++)
                printf("%d ", a[i][j]);
            printf("%d\n", a[i][sz - 1]);
        }
        printf("\n");
    }
};

Matrix pow_mod(Matrix a, int p) {
    Matrix ans;
    for(int i = 0; i < sz; i++) ans.a[i][i] = 1;
    while(p) {
        if(p & 1) ans = ans * a;
        a = a * a;
        p >>= 1;
    }
    return ans;
}

int main()
{
    while(scanf("%d%d", &n, &k) == 2 && n) {
        Matrix M;
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) {
                scanf("%d", &a[i][j]);
                a[i][j] %= MOD;
                M.a[i][j] = a[i][j];
            }

        if(!k) {
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n - 1; j++)
                    printf("%d ", i == j ? 1 : 0);
                printf("%d\n", i == n - 1 ? 1 : 0);
            }
            printf("\n");
            continue;
        }
        
        for(int i = n; i < n * 2; i++)
            M.a[i][i] = M.a[i][i - n] = 1;

        sz = n * 2;
        M = pow_mod(M, k);
        memset(b, 0, sizeof(b));
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) if(M.a[i + n][j])
                for(int k = 0; k < n; k++)
                    add(b[i][k], mul(M.a[i + n][j], a[j][k]));

        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n - 1; j++)
                printf("%d ", b[i][j]);
            printf("%d\n", b[i][n - 1]);
        }
        printf("\n");
    }

    return 0;
}

方法二

有如下递归式：

\(S_k=(E+A^{\frac{k}{2}})S^{\frac{k}{2}}\)，k是偶数
\(S_k=(E+A^{\frac{k}{2}})S^{\frac{k}{2}}+A^k\)，k是奇数

所以也可以直接递归求解。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int maxn = 40;
const int MOD = 10;

int mul(int a, int b) { return a * b % MOD; }

void add(int& a, int b) { a += b; if(a >= MOD) a -= MOD; }

int n, k;

struct Matrix
{
    int a[maxn][maxn];

    Matrix() { memset(a, 0, sizeof(a)); }
    
    void E() {
        memset(a, 0, sizeof(a));
        for(int i = 0; i < maxn; i++) a[i][i] = 1;
    }
    
    Matrix operator + (const Matrix& t) const {
        Matrix ans;
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) {
                add(ans.a[i][j], a[i][j]);
                add(ans.a[i][j], t.a[i][j]);
            }
        return ans;
    }
    
    Matrix operator * (const Matrix& t) const {
        Matrix ans;
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) if(a[i][j])
                for(int k = 0; k < n; k++)
                    add(ans.a[i][k], mul(a[i][j], t.a[j][k]));
        return ans;
    }

    void output() {
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n - 1; j++)
                printf("%d ", a[i][j]);
            printf("%d\n", a[i][n - 1]);
        }
        printf("\n");
    }
};

Matrix pow_mod(Matrix a, int p) {
    Matrix ans;
    ans.E();
    while(p) {
        if(p & 1) ans = ans * a;
        a = a * a;
        p >>= 1;
    }
    return ans;
}

Matrix E;

Matrix sum(Matrix a, int p) {
    if(p == 1) return a;
    Matrix ans;
    ans = (E + pow_mod(a, p / 2)) * sum(a, p / 2);
    if(p & 1) ans = ans + pow_mod(a, p);
    return ans;
}

int main()
{
    E.E();
    while(scanf("%d%d", &n, &k) == 2) {
        if(n == 0 && k == 0) break;
        Matrix a;
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++) {
                scanf("%d", &a.a[i][j]);
                a.a[i][j] %= 10;
            }

        if(k == 0) {
            E.output();
            continue;
        }

        a = sum(a, k);
        a.output();
    }

    return 0;
}

转载于:https://www.cnblogs.com/AOQNRMGYXLMV/p/5260771.html