[WOJ26 Lost in WHU]矩阵快速幂

分类：Matrix Math

1. 题目链接

[WOJ26 Lost in WHU]

2. 题意描述

给定一个 N 个顶点的对称矩阵G (1≤N≤100) ， Aij 表示可以从点 i 一步到点j，求在 T 时刻之前能从1到 N 的方案数。(1≤T≤109) 注意：如果在时刻 tx 已经到达了点 N ，则终止。

3. 解题思路

题意是求在0,1,2,3,…,T时刻到达 N 点的方案之和。 可以这么考虑，求经过tx步达到的 N 点的方案数可以这么计算，tx−1步是在 1,2,3,…,N−1 这几个顶点中走，最后一步再走到 N 点。 那么，问题就很简单的转化为熟悉的矩阵快速幂了，经过tx步达到的 N 点的方案数就是除去点N的邻接矩阵的 tx−1 次方，然后加上最后一步的情况。 问题现在转化为对矩阵幂求和了，即求矩阵 ∑T−1i=0Ai 。 设 Sx=∑x−1i=0Ai 。 下面是求 Sx 的3种做法：

分治法： 利用 Sx=Sx2+Sx2∗Ax2 ，进行分治即可。复杂度： O(N3log2(T)) 构造矩阵法：根据 Sx=A0+Sx−1∗A 很容易构造出矩阵转移式： [SxE]=[AOEE]∗[Sx−1E] E 是单位矩阵，O是零矩阵。这样搞的复杂度是 O((2N)3log(T)) 通项公式法?： 这个做法是我yy的。理论上应该是可行的。类比等比数列求和。 Sx=E−AxE−A 。这样就需要求一个矩阵的逆（高斯消元来做，好麻烦）。

最后这个题目矩阵挺大的。递归一多肯定爆栈，所以可以必须用全局变量或静态变量。推荐用静态变量，返回引用值来做，简单易实现，效率还高。

4. 实现代码

/** 法一：构造矩阵法 **/ #include <bits/stdc++.h> using namespace std; typedef long long LL; typedef long double LB; typedef pair<int, int> PII; typedef pair<LL, LL> PLL; typedef vector<int> VI; const int INF = 0x3f3f3f3f; const LL INFL = 0x3f3f3f3f3f3f3f3fLL; const double eps = 1e-8; const double PI = acos(-1.0); const int MAXN = 100000 + 5; const int MX = 100 + 1; const int MOD = 1e9 + 7; int n, m, t, msz; struct Mat { LL d[MX << 1][MX << 1]; void I() { memset(d, 0, sizeof(d)); } void E() { I(); for(int i = 1; i <= msz; i++) d[i][i] = 1; } Mat& operator * (const Mat& e) const { static Mat ret; ret.I(); for(int k = 1; k <= msz; k++) { for(int i = 1; i <= msz; i++) { if(d[i][k] == 0) continue; for(int j = 1; j <= msz; j++) { ret.d[i][j] += d[i][k] * e.d[k][j]; ret.d[i][j] %= MOD; } } } return ret; } Mat& operator ^ (int b) const { static Mat x(*this), ret; ret.E(); while(b > 0) { if(b & 1) ret = ret * x; x = x * x; b >>= 1; } return ret; } void P() const { for(int i = 1; i <= msz; i++) { for(int j = 1; j <= msz; j++) { printf("%lld ", d[i][j]); } puts(""); } } } init, tran, ans; int G[MX]; int main() { #ifdef ___LOCAL_WONZY___ freopen("input.txt", "r", stdin); #endif // ___LOCAL_WONZY___ static int u, v; scanf("%d %d", &n, &m); memset(G, 0, sizeof(G)); init.I(); tran.I(); msz = ((n - 1) << 1); for(int i = 1; i <= n - 1; ++i) { init.d[i][i] = 1; tran.d[i][i + n - 1] = 1; } for(int i = n; i <= msz; ++i) { init.d[i][i - n + 1] = 1; tran.d[i][i] = 1; } for(int i = 1; i <= m; ++i) { scanf("%d %d", &u, &v); if(u > v) swap(u, v); if(v == n) { G[u] = 1; continue; } tran.d[u][v] = ++ tran.d[v][u]; } scanf("%d", &t); tran = (tran ^ (t - 1)); ans = tran * init; LL rs = 0; for(int i = 1; i <= (n - 1); ++i) { if(!G[i]) continue; rs = (rs + ans.d[1][i] * G[i]) % MOD; } printf("%lld\n", rs); #ifdef ___LOCAL_WONZY___ cout << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC * 1000 << " ms." << endl; #endif // ___LOCAL_WONZY___ return 0; } /** 法二：分治法 **/ #include <bits/stdc++.h> using namespace std; typedef long long LL; typedef long double LB; typedef pair<int, int> PII; typedef pair<LL, LL> PLL; typedef vector<int> VI; const int INF = 0x3f3f3f3f; const LL INFL = 0x3f3f3f3f3f3f3f3fLL; const double eps = 1e-8; const double PI = acos(-1.0); const int MAXN = 100000 + 5; const int MX = 100 + 1; const int MOD = 1e9 + 7; int n, m, t; struct Mat { LL d[MX][MX]; void I() { memset(d, 0, sizeof(d)); } void E() { I(); for(int i = 1; i <= n - 1; i++) d[i][i] = 1; } Mat& operator * (const Mat& e) const { static Mat ret; ret.I(); for(int k = 1; k <= n - 1; k++) { for(int i = 1; i <= n - 1; i++) { if(d[i][k] == 0) continue; for(int j = 1; j <= n - 1; j++) { ret.d[i][j] += d[i][k] * e.d[k][j]; ret.d[i][j] %= MOD; } } } return ret; } Mat& operator ^ (int b) const { static Mat x, ret; x = (*this), ret.E(); while(b > 0) { if(b & 1) ret = ret * x; x = x * x; b >>= 1; } return ret; } Mat& operator + (const Mat& e) const { static Mat ret; for(int i = 1; i <= n - 1; i++) { for(int j = 1; j <= n - 1; j++) { ret.d[i][j] = (d[i][j] + e.d[i][j]) % MOD; } } return ret; } void P() const { for(int i = 1; i <= n - 1; i++) { for(int j = 1; j <= n - 1; j++) { printf("%lld ", d[i][j]); } puts(""); } } } init, tran, ans; int G[MX]; Mat& ask(const Mat& a, int b) { static Mat ret; if(b <= 1) { ret.E(); return ret; } int md = (b >> 1); ret = ask(a, md); ret = ret + ret * (a ^ md); if(b & 1) ret = ret + (a ^ (b - 1)); return ret; } int main() { #ifdef ___LOCAL_WONZY___ freopen("input.txt", "r", stdin); #endif // ___LOCAL_WONZY___ static int u, v; scanf("%d %d", &n, &m); memset(G, 0, sizeof(G)); tran.I(); for(int i = 1; i <= m; ++i) { scanf("%d %d", &u, &v); if(u > v) swap(u, v); if(v == n) { G[u] = 1; continue; } tran.d[u][v] = ++ tran.d[v][u]; } scanf("%d", &t); ans = ask(tran, t); LL rs = 0; for(int i = 1; i <= (n - 1); ++i) { if(!G[i]) continue; rs = (rs + ans.d[1][i] * G[i]) % MOD; } printf("%lld\n", rs); #ifdef ___LOCAL_WONZY___ cout << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC * 1000 << " ms." << endl; #endif // ___LOCAL_WONZY___ return 0; }