大禹治水题解

出题

---

当 $Case\in[1,10]$ 时，由于数据范围很小，完全可以暴力深搜计算到达每个节点的方案数，最后将与节点 $1$ 相邻的方案数相加取余就是答案，时间复杂度为 $O(n^n)$，这种做法可以通过这部分数据。

---

对于统计图中路径方案数，可以先创建一个邻接矩阵 M[i][j]（下标从 $1$ 开始），M[i][j]=M[j][i]=x 表示节点 $i$ 与节点 $j$ 之间有 $x$ 条无向边，注意当 $u=v$ 时，M[u][v] 只能被计算一次；由于大禹不能走回家里，只能有以 $1$ 为起点的边，不能有以 $1$ 为终点的边，所以 G[i][1] 需要手动清零。

两个邻接矩阵相乘，得到的第 $i$ 行第 $j$ 列的值，就是从节点 $i$ 出发，走两步能到达节点 $j$ 的方案数，题目要求 $k$ 步以内的方案数总和，即求以下公式的结果：

$$M+M^2+M^3+\cdots M^k$$
取结果矩阵的第 $1$ 行所有最初与节点 $1$ 相邻的节点位置的方案数总和，即为答案。

---

当 $Case\in[11,30]$ 时，由于 $k$ 非常大，每次矩阵乘法的时间复杂度都为 $O(n^3)$，若使用顺序乘法，总共需要计算 $50^3\times10^6=1.25\times10^{11}$ 次，必然超时，需要优化。

考虑计算 $sum_k=\sum_{i=1}^k M^i$，有以下递推式：

$$sum_k=M\times(sum_{k-1}+I)$$
其中 $I$ 为单位矩阵，$M$ 为邻接矩阵，将其转化为矩阵运算（矩阵内每一项都是一个 $n^2$ 的矩阵）：
$$\begin{pmatrix}sum_k\\M\end{pmatrix}=\begin{pmatrix}M&I\\0&I\end{pmatrix}\begin{pmatrix}sum_{k-1}\\M\end{pmatrix}=\begin{pmatrix}M&I\\0&I\end{pmatrix}^k\begin{pmatrix}0\\M\end{pmatrix}$$
由于 $k$ 次幂的每一项都相同，且矩阵乘法满足结合律，因此可以使用矩阵快速幂在 $O(\log_2 k)$ 次内计算出矩阵 $sum_k$，该解法的时间复杂度为 $O(2^3n^3\log_2 k)$。

---

当 $Case\in[31,50]$ 时，$k$ 的数据范围达到 $10^{100000}$，上一种解法已不适用，观察到节点 $n$ 与膜数 $p$ 的值非常小，考虑计算循环节，由于矩阵 $M$ 第一列元素全为 $0$，之后无论如何相乘第一列都必然为 $0$，因此可以忽略第一列的状态变化。

除去第 $1$ 列共有 $n(n-1)$ 个元素在不断变化，可知循环节长度必然小于等于 $p^{n(n-1)}$ ，因此最多只要计算出前 $p^{n(n-1)}$ 项 $sum_k$ 即可。

循环节状态必然为：

$$sum_0\rightarrow sum_1\rightarrow sum_2\rightarrow\cdots\rightarrow\underline{sum_a}\rightarrow sum_{a+1}\rightarrow\cdots\rightarrow sum_b\rightarrow\underline{sum_a}$$
此时循环节长度为 $len=b+1-a$，若 $k<a$，则直接对应 $sum_k$ 的结果，若 $k\geq a$，则 $sum_k$ 将对应到 $sum_{(k-a)\%len+a}$ 的结果。

由于 $k$ 非常大，但 $k$ 对 $len$ 取余前后再对 $a$ 做减法的结果是一致的，因此可以不必做大数减法计算，先用 $k$ 对 $len$ 取余后再对 $a$ 做减法同样可以得到正确的循环节下标。

时间复杂度为 $O(n^3\times p^{n(n-1)}+\log_{10}k)$，最坏情况下需要执行约 $3.4\times10^7$ 次计算。

暴力解法

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn = 100;
int T, n, m, k, p, ans, u, v;
int cnt[maxn];
bool isDoorway[maxn];
vector<int> G[maxn];
void dfs(int root, int k) {
    ++cnt[root];
    if (k == 0) {
        return ;
    }
    for (auto pos: G[root]) {
        if (pos == 1) {
            continue;
        }
        dfs(pos, k - 1);
    }
}
int main() {
    cin >> T;
    cin >> n >> m >> k >> p;
    for (int i = 0; i < m; ++i) {
        cin >> u >> v;
        if (u > v) {
            swap(u, v);
        }
        G[u].push_back(v);
        if (u == 1) {
            isDoorway[v] = true;
        } else {
            if (u != v) {
                G[v].push_back(u);
            }
        }
    }
    dfs(1, k);
    for (int i = 2; i <= n; ++i) {
        if (isDoorway[i]) {
            ans += cnt[i];
        }
    }
    cout << ans % p << endl;
    return 0;
}

标程

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int matrixSize = 50;
const int maxn = 100000 + 100;
int MOD;
bool isDoorway[50];
LL add(LL a, LL b) {
    a += b;
    if(a >= MOD) {
        return a - MOD;
    }
    return a;
}
struct Matrix {
    int size;
    LL num[matrixSize][matrixSize];
    const Matrix& operator=(const Matrix &m) {
        size = m.size;
        for(int i = 0; i < size; ++i) {
            memcpy(num[i], m.num[i], sizeof(LL) * size);
        }
        return (*this);
    }
    void init() {
        for(int i = 0; i < size; ++i) {
            memset(num[i], 0, sizeof(LL) * size);
            num[i][i] = 1;
        }
    }
    void clear() {
        for (int i = 0; i < size; ++i) {
            memset(num[i], 0, sizeof(LL) * size);
        }
    }
    void set(int n, int m) {
        size = n;
        int u, v;
        for (int i = 0; i < m; ++i) {
            scanf("%d%d", &u, &v);
            --u;
            --v;
            ++num[u][v];
            if (u != v) {
                ++num[v][u];
            }
        }
        for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
                if (j == 0) {
                    num[i][j] = 0;
                } else {
                    if (i == 0) {
                        isDoorway[j] = (num[i][j] != 0);
                    }
                    num[i][j] %= MOD;
                }
            }
        }
    }
    const Matrix& operator*(const Matrix &m) const {
        static Matrix ans;
        ans.size = m.size;
        for(int i = 0; i < size; ++i) {
            for(int j = 0; j < size; ++j) {
                ans.num[i][j] = 0;
                for(int k = 0; k < size; ++k) {
                    ans.num[i][j] = add(ans.num[i][j], num[i][k] * m.num[k][j] % MOD);
                }
            }
        }
        return ans;
    }
    const Matrix& operator*=(const Matrix &m) {
        (*this) = (*this) * m;
        return *this;
    }
    const Matrix& operator+(const Matrix &m) const {
        static Matrix ans;
        ans.size = m.size;
        for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
                ans.num[i][j] = add(num[i][j], m.num[i][j]);
            }
        }
        return ans;
    }
    const Matrix& operator+=(const Matrix &m) {
        for (int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
                num[i][j] = add(num[i][j], m.num[i][j]);
            }
        }
        return *this;
    }
    LL hash() {
        LL ret = 0;
        for (int i = 0; i < size; ++i) {
            for (int j = 1; j < size; ++j) {
                ret = ret * 10 + num[i][j];
            }
        }
        return ret;
    }
};
Matrix one, zero;
struct MMatrix {
    int size;
    Matrix m[2][2];
    void set(const Matrix &m) {
        size = 2;
        this->m[0][0] = m;
        this->m[0][1] = one;
        this->m[1][0] = zero;
        this->m[1][1] = one;
    }
    const MMatrix& operator=(const MMatrix &mm) {
        size = mm.size;
        for(int i = 0; i < size; ++i) {
            for (int j = 0; j < size; ++j) {
                m[i][j] = mm.m[i][j];
            }
        }
        return (*this);
    }
    void init() {
        for(int i = 0; i < size; ++i) {
            m[i][i] = one;
            m[i][i ^ 1] = zero;
        }
    }
    const MMatrix& operator*=(const MMatrix &mm) {
        static MMatrix ans;
        ans.size = size;
        for(int i = 0; i < size; ++i) {
            for(int j = 0; j < size; ++j) {
                ans.m[i][j] = zero;
                for(int k = 0; k < size; ++k) {
                    ans.m[i][j] += m[i][k] * mm.m[k][j];
                }
            }
        }
        (*this) = ans;
        return ans;
    }
    void fastPow(LL n) {
        static MMatrix ans;
        ans.size = size;
        for(ans.init(); n != 0; n >>= 1) {
            if((n & 1) == 1) {
                ans *= (*this);
            }
            (*this) *= (*this);
        }
        (*this) = ans;
    }
};
int T, n, m, k, ans;
Matrix matrix, matrixSum;
MMatrix mmatrix;
char kStr[maxn];
LL h;
LL pow10[20];
map<LL, int> hashIndex;
vector<LL> hs;
int getNum(LL h, int x, int y) {
    int r = n - 1 - x;
    int c = n - 1 - y;
    return h / pow10[r * (n - 1) + c] % 10;
}
void solve1() {
    scanf("%d%d%d%d", &n, &m, &k, &MOD);
    matrix.set(n, m);
    one.size = n;
    one.init();
    zero.size = n;
    zero.clear();
    mmatrix.set(matrix);
    mmatrix.fastPow(k);
    mmatrix.m[0][1] *= matrix;
    ans = 0;
    for (int i = 1; i < n; ++i) {
        if (isDoorway[i]) {
            ans = (ans + mmatrix.m[0][1].num[0][i]) % MOD;
        }
    }
    printf("%d\n", ans);
}
void solve2() {
    pow10[0] = 1;
    for (int i = 1; i < 18; ++i) {
        pow10[i] = pow10[i - 1] * 10;
    }
    scanf("%d%d%s%d", &n, &m, kStr, &MOD);
    matrix.set(n, m);
    matrixSum.size = n;
    matrixSum.clear();
    h = matrixSum.hash();
    do {
        hashIndex[h] = hs.size();
        hs.push_back(h);
        matrixSum = matrixSum * matrix + matrix;
        h = matrixSum.hash();
    } while (hashIndex.find(h) == hashIndex.end());
    int loopIndex = hashIndex[h];
    int loopLength = hs.size() - loopIndex;
    int index = 0;
    for (int i = 0; kStr[i] != '\0'; ++i) {
        index = (index * 10 + (kStr[i] - '0')) % loopLength;
    }
    index = ((index - loopIndex) % loopLength + loopLength) % loopLength + loopIndex;
    for (int j = 1; j < n; ++j) {
        if (isDoorway[j]) {
            ans = (ans + getNum(hs[index], 0, j)) % MOD;
        }
    }
    printf("%d\n", ans);
}
int main() {
    scanf("%d", &T);
    if (T <= 30) {
        solve1();
        return 0;
    }
    solve2();
    return 0;
}

数据生成

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int testCases = 50;
ofstream fout;
int T, n, m, k, p, u, v;
string numToStr(int x) {
    string str;
    while (x != 0) {
        str = char('0' + x % 10) + str;
        x /= 10;
    }
    return str;
}
string getFileName(int index) {
    return "data/" + numToStr(index) + ".in";
}
int getRand(int l, int r) {
    return (LL)rand() * rand() * rand() % (r - l + 1) + l;
}
void writeToFile() {
    for (int i = 0; i < m; ++i) {
        u = getRand(1, n);
        v = getRand(1, n);
        fout << u << " " << v << endl;
    }
}
int main() {
    srand(time(0));
    // data in 20%, 2 <= n,m,k <= 7, p = 1e9+7
    int one = testCases * 20 / 100;
    for (int i = 1; i <= one; ++i) {
        fout.open(getFileName(i));
        fout << i << endl;
        n = getRand(5, 7);
        m = getRand(5, 7);
        k = getRand(5, 7);
        p = 1000000007;
        fout << n << " " << m << " " << k << " " << p << endl;
        writeToFile();
        fout.close();
    }
    // data in 60%, 2 <= n <= 50, 1 <= m <= 1e6, 1 <= k <= 1e6, 1e7 <= p <= 1e9
    int two = testCases * 60 / 100;
    for (int i = one + 1; i <= two; ++i) {
        fout.open(getFileName(i));
        fout << i << endl;
        n = getRand(45, 50);
        m = getRand(100000, 1000000);
        k = getRand(999900, 1000000);
        p = getRand(10000000, 1000000000);
        fout << n << " " << m << " " << k << " " << p << endl;
        writeToFile();
        fout.close();
    }
    // data in 100%, 2 <= n <= 4, 1 <= m <= 1e6, 1 <= k <= 1e100000, 1 <= p <= 3
    for (int i = two + 1; i <= testCases; ++i) {
        fout.open(getFileName(i));
        fout << i << endl;
        n = getRand(3, 4);
        m = getRand(100000, 1000000);
        p = getRand(2, 3);
        int lenOfk = getRand(50000, 100000);
        fout << n << " " << m << " ";
        fout << getRand(1, 9);
        for (int j = 1; j < lenOfk; ++j) {
            fout << getRand(0, 9);
        }
        fout << " " << p << endl;
        writeToFile();
        fout.close();
    }
    return 0;
}