矩阵快速幂-HDU-07-1008 循环图

题意

坎格鲁斯普雷被困在了一张循环图里，这张循环图有无数个节点，初始时坎格鲁斯普雷在 号节点。

循环图是边存在着一定循环关系的图，循环图里面的边可以用循环周期 和 对三元组 表示。每对三元组 表示，对于循环图内所有的满足 的点对 ，都存在有 ​ 条从点 通往点 的边。

现在，坎格鲁斯普雷知道了这张循环图的第 个节点到第 个节点各存在着一个出口，坎格鲁斯普雷需要到达这些节点中的任意一个才能逃出循环图（到达有出口存在的节点后不一定要立刻逃出）。坎格鲁斯普雷想请你帮他算算，他有多少种逃出这张循环图的方式。由于答案可能很大，你需要输出答案对 取模后的结果。

朴素思路

首先可以想到朴素 ，表示从走到有种方案，转移可以 = 。但是这样显然会超时。

矩阵快速幂思路

由于、很大，很小，而又存在非常显然的递推关系，即，十分符合矩阵的乘法，可以想到矩阵快速幂。
现在就是要考虑矩阵需要记录什么。由于需要求区间，可以转化成前缀和的差，所以矩阵需要记录前缀。而前面提到的递推关系也很自然的想到用矩阵表示从走到的方案数。

然后就是怎么递推了。我们发现 ，的范围很好，但是会到的范围，我们不妨先用矩阵记录的方案数，再用矩阵记录从 一步走到的方案数，这样利用矩阵乘法就可以算出从 走到的总方案数（总方案数意思是可以有中转，比如1->2->10）。
再考虑前缀和。这肯定是要再开一维的。考虑到要利用到矩阵快速幂，我们发现两个极好的性质：
1. 由题意，这样的点对之间有边，所以既可以表示以为起点，走到的方案数，又可以表示从开始，走到的方案数。你会发现就是从开始，走到的方案数！
2. 由于前缀和要和方案数存在同一个矩阵里面，那么如果前缀和与方案数在同一个维度（意思是说，矩阵表示从走到，前缀和也表示走到的方案数，这两个都是(k + 1)维度的），那么递推肯定是不好推的，因为我现在还没有求得下一维度的方案数，我根本就没法得到下一个维度的前缀和。说到这里就很有感触了，那就是把前缀和降一个维度！这样一个矩阵里面存的是这个维度的方案数和上一个维度的前缀和，而上一维度的前缀和+这一维度的方案是就是这一维度的前缀和，也就是下一矩阵中的前缀和！

注意矩阵的初始化就AC了（

code

#include<bits/stdc++.h>
using namespace std;
#define endl '\n'
using ll = long long;
const int mod = 1e9 + 7;

struct matrix {
    vector<vector<ll>> t;
    int n, m;

    matrix(int _n, int _m) {
        n = _n, m = _m;
        t.resize(n + 1);
        for (int i = 1; i <= n; i++) {
            t[i].resize(m + 1);
        }
    }

    matrix operator *(const matrix &a) const {
        matrix res(n, a.m);
        for (int i = 1; i <= res.n; i++) {
            for (int j = 1; j <= res.m; j++) {
                for (int k = 1; k <= m; k++) {
                    res.t[i][j] += t[i][k] * a.t[k][j];
                    res.t[i][j] %= mod;
                }
            }
        }
        return res;
    }

    matrix m_qpow(ll k) {
        matrix res(n, n);
        auto A = *this;
        for (int i = 1; i <= n; i++) {
            res.t[i][i] = 1;
        }
        while (k) {
            if (k & 1) {
                res = res * A;
            }
            A = A * A;
            k >>= 1;
        }
        return res;
    }
};

void solve() {
    int n, m;
    ll L, R;
    cin >> n >> m >> L >> R;
    vector<vector<ll>> ed(n + 1, vector<ll>(n + 1)), W(n + 1, vector<ll>(n + 1)), dp(n + 1, vector<ll>(n + 1));
    vector<bool> vis(n + 1);
    for (int i = 1; i <= m; i++) {
        int u, v, w;
        cin >> u >> v >> w;
        if (v <= n) {
            ed[u][v] = w;
        }
        else {
            W[u][v - n] = w;
        }
        // 分别记录，这样好处理出D和Q矩阵
    }

    auto dfs = [&](auto self, int u) -> void {
        if (vis[u]) return;
        vis[u] = true;
        dp[u][u] = 1;
        for (int v = u + 1; v <= n; v++) {
            for (int k = u + 1; k <= v; k++) {
                self(self, k);
                dp[u][v] = (dp[u][v] + 1ll * ed[u][k] * dp[k][v] % mod) % mod;
            } 
        }
    };

    for (int i = 1; i <= n; i++) {
        dfs(dfs, i);
        // 预处理D矩阵
    }

    auto get = [&](ll x) -> ll {
        ll k = x / n;
        // 这里要注意因为矩阵是n个一推的，所以要留下x % n个直接加
        matrix D(n + 1, n + 1);
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                D.t[i][j] = dp[i][j];
            }
        }

        D.t[n + 1][n + 1] = 1;

        matrix Q(n + 1, n + 1);
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                Q.t[i][j] = W[i][j];
            }
        }
        for (int i = 1; i <= n + 1; i++) {
            Q.t[i][n + 1] = 1;
            // 给第n + 1维赋值，好直接处理出A的前缀和列（第n + 1列）
        }

        auto A = D * Q;
        auto res = A.m_qpow(k);
        ll ans = res.t[1][n + 1];
        res = res * D;
        // *A实际上是*D*Q，即先自己内部计算然后再递推到下一个，这里算剩下的不需要递推到下一个，所以只要*D即可
        for (int i = 1; i <= x % n; i++) {
            ans = (ans + res.t[1][i]) % mod;
        }
        return ans;
    };

    ll ans = get(R) - get(L - 1);
    // 前缀和相减得到区间和
    ans = (ans + mod) % mod;
    cout << ans << endl;
}


int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout.tie(0);
    int tt = 1;
    cin >> tt;
    while (tt--) solve();
    return 0;
}