题解 UVA11149 Power of Matrix

简要题意：

给定一个 $n \times n$ 的矩阵 $A$ 和一个数 $k$，求 $\sum_{i = 1}^{k}A^i $。

法一：暴力求解，单次时间复杂度 $O(n^3k)$

法二：分治，单次时间复杂度 $O(n^3log^2k)$

我们可以按照分治的定义照猫画虎，将原式变成这样：
$$A+A^2+A^3+A^4+\ldots+A^K = (A+A^2+A^3+\ldots+A^{\frac{k}{2}})+A^{\frac{k}{2}}(A+A^2+A^3+\ldots+A^{\frac{k}{2}})$$

然后递归处理即可，分治部分的代码大概长这样：

1
2
3
4
5
6
7
8
mat solve(mat bas, LL pos) { \\bas是原矩阵，pos代表当前k的大小
    if (pos == 1) { return bas; }
    mat s; clear(s);
    for (int i = 0; i < n; ++i) { s.m[i][i] = 1; }
    s = (s + qpow(bas, pos >> 1)) * solve(bas, pos >> 1);
    if (pos & 1) { s = s + qpow(bas, pos); }
    return s;
}

于是我们惊喜的发现，这种做法又慢又不稳定，所以我们考虑法三。

法三：倍增，单次时间复杂度 $O(n^3logk)$

我们观察要求的这个式子，发现 $i$ 要枚举到 $k$ 非常慢，那么我们就想办法优化这个枚举的过程，也就是对 $i$ 进行倍增，则有：
$$\sum_{i=0}^{log_2k} A^{2^i}(\sum_{j=1}^{2^i} A^j)$$
我们显然不能每次都从头到尾加一遍，所以我们用两个数组 $f$ 和 $s$对这个式子进行递推。
$$f_i=A^{2^i},s_i=\sum_{j=1}^{2^i}A^j$$

不难发现在递推 $f_i$ 时只要每次把 $f_{i-1}$ 平方就可以了，而 $s_i$ 略微复杂一些，求它其实很像刚刚分治的做法，即 $s_i=s_{i-1}+f_{i-1} \times s_{i-1}$。

如果题目高速我们 $k$ 保证是 2 的若干次幂，那么这个题已经做完了，但很遗憾我们没有这个条件，所以还要想办法处理剩余的部分。本过程较为复杂，建议多看几遍。 我们知道整数的一个重要性质就是可以将其拆解为若干 2 的幂之和，假设我们的 $i$ 最后枚举到了 $m$ 使得 $2^m > k$ 且 $2^{m-1} < k$，那么此时我们就想要拆解 $k-2^{m-1}$ 这一部分，那么显然我们最简单的方法是从 $m-1$ 开始向下枚举，设循环变量为 $i$，当枚举到一个 $i$ 使得 $2^{m-1}+2^i<k$ 的时候，就从刚刚处理出来的两个数组中取出第 $i$ 项加到答案里就可以了。这样枚举到 $s^{m-1}+2^i=k$ 时说明我们已经完成对答案的处理，输出答案就可以了。

下面是参考代码（倍增过程在 $\operatorname{solve}$ 函数部分）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#define clear(q) memset(q.m, 0, sizeof(q.m))
#define LL long long
using namespace std;

namespace SHAWN {
    LL n, k;
    const int N = 41, M = 64, mod = 10;
    struct mat { LL m[N][N]; };
    mat f[M], s[M];
    mat operator* (const mat &x, const mat &y) {
        mat ans;
        clear(ans);
        for (int i = 0; i < n; ++i) {
            for (int k = 0; k < n; ++k) {
                if (x.m[i][k]) {
                    for (int j = 0; j < n; ++j) {
                        ans.m[i][j] = (ans.m[i][j] % mod + (x.m[i][k] * y.m[k][j]) % mod) % mod;
                    }
                }
            }
        }
        return ans;
    }
    mat operator+ (const mat &x, const mat &y) {
        mat ans;
        clear(ans);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                ans.m[i][j] = (x.m[i][j] + y.m[i][j]) % mod;
            }
        }
        return ans;
    }
    mat qpow(mat a, LL b) {
        mat z;
        clear(z);
        for (int i = 0; i < n; ++i) { z.m[i][i] = 1; }
        for (int i = 0; i < n; ++i) {
            if (b & 1) { z = z * a; }
            a = a * a; b >>= 1;
        }
        return z;
    }
    void output(mat o) {
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                cout << o.m[i][j] % mod;
                if (j != n - 1) { cout << ' '; }
            }
            cout << '\n'; 
        } 
        cout << '\n';
    }
    void solve() {
        LL pos = 2, p = 1; // pos用来枚举2的次幂，p用来枚举2的幂次，pos的指数即为p，二者同加同减
        while (pos <= k) {
            f[p] = f[p - 1] * f[p - 1];
            s[p] = s[p - 1] + s[p - 1] * f[p - 1];
            pos <<= 1; ++p;
        }
        pos >>= 1; --p; // 第一次枚举完后肯定会超出k的范围，所以取上一次枚举到的幂
        LL posback = pos; // posback表示拆分剩余的数时需要用到的2的次幂
        mat po = f[p], ans = s[p]; // po表示当前已经枚举到的幂次，ans为答案，初始值为第一次枚举时枚举到最大的值
        for (int i = p; i >= 0; --i) {
            if (pos + posback <= k) { 
                ans = ans + po * s[i]; // 加到答案里
                posback += pos;
                po = po * f[i];
            }
            pos >>= 1; // pos的幂次要和p同步减一
        }
        output(ans); // 处理完输出
    }
    int work()
    {
        while (cin >> n >> k && n && k) {
            for (int i = 0; i < n; ++i) {
                for (int j = 0; j < n; ++j) {
                    cin >> s[0].m[i][j];
                    f[0].m[i][j] = s[0].m[i][j] = s[0].m[i][j] % mod; 
                }
            }
            solve();
        }
        return 0;
    }
}

signed int main() {
    ios :: sync_with_stdio(false);
    cin.tie(nullptr);
    return SHAWN :: work();
}

ps. 特别提醒，UVa要注意输入和输出格式

说些什么吧！

twikoo