WOJ26 矩阵快速幂

WOJ26 矩阵快速幂：给定一个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种做法：

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

4. 实现代码

/** 法一：构造矩阵法 **/
#include 
using namespace std;

typedef long long LL;
typedef long double LB;
typedef pair PII;
typedef pair PLL;
typedef vector 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 
using namespace std;

typedef long long LL;
typedef long double LB;
typedef pair PII;
typedef pair PLL;
typedef vector 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;
}