[WOJ26 Lost in WHU]矩阵快速幂

来源：互联网发布：leqee网络上海有公司吗编辑：程序博客网时间：2024/06/01 19:38

[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很容易构造出矩阵转移式： $[S x E] = [A O E E] * [S x - 1 E]$ 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;}

1 0