[WOJ26 Lost in WHU]矩阵快速幂

[WOJ26 Lost in WHU]矩阵快速幂

分类：Matrix Math

1. 题目链接

[WOJ26 Lost in WHU]

2. 题意描述

给定一个$N$个顶点的对称矩阵$G$$(1\le N\le 100)$，$A_{ij}$表示可以从点$i$一步到点$j$，求在$T$时刻之前能从$1$到$N$的方案数。$(1\le T \le 10^9)$
注意：如果在时刻$t_x$已经到达了点$N$，则终止。

3. 解题思路

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

• 分治法： 利用$S_x=S_{\frac{x}{2}}+S_{\frac{x}{2}}*A^{\frac{x}{2}}$，进行分治即可。复杂度：$O(N^3log^2(T))$
• 构造矩阵法：根据$S_x=A^0+S_{x-1}*A$很容易构造出矩阵转移式： $$\begin{bmatrix}S_x\\E\end{bmatrix}=\begin{bmatrix}A & E\\O &E\end{bmatrix}*\begin{bmatrix}S_{x-1}\\E\end{bmatrix}$$ $E$是单位矩阵，$O$是零矩阵。这样搞的复杂度是$O((2N)^3log(T))$
• 通项公式法?： 这个做法是我yy的。理论上应该是可行的。类比等比数列求和。$S_x=\frac{E-A^x}{E-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;
}