跳转至

状压 DP

枚举子集

有时，对于一个二进制状态 \(i\)，我们只希望从它的子集转移，因为这样能有效的减少时间复杂度。

Y555 最优组队

给定 \(n\) 个人（\(n\le 16\)），你需要把他们分成若干组，满足所有组的权值和最大。题目给定所有可能出现的组对应的权值。

我们注意到 \(n\) 很小，可以指数复杂度通过，因此考虑状压。设二进制状态 \(i\)，表示对应集合中的人都已组好队，最大的权值和。很容易写出状态转移方程：

\[ f_i=\max_{j\subset i}\{f_j+w_{i-j}\} \]

对于状态 \(i\)，暴力枚举所有状态 \(j\)，再判断是否满足 \(j\subset i\)，决定是否转移。时间复杂度为 \(O(4^n)\)，“不能”通过本题。

1
2
3
4
5
6
7
8
9
for(int i = 1; i <= (1 << n) - 1; i++) {
    for(int j = 0; j < i; j++) {
        if((i | j) == i){
            f[i] = max(f[i], f[j] + a[i ^ j]);
        }
    }
}

cout << f[(1 << n) - 1] << endl;

我也不知道这段暴力为什么就水过了，可能是带一个小常数。

我们注意到，大量的时间被浪费在了枚举那些不是 \(i\) 的子集的 \(j\) 上面。在此我们提供一种可以只遍历到所有真子集的方法：

1
2
3
4
5
6
7
for(int i = 1; i <= (1 << n) - 1; i++) {
    int j = i;
    while(j) {
        j = (j - 1) & i;
        f[i] = max(f[i], f[j] + a[i ^ j]);
    }
}

通过此种优化，时间复杂度降为 \(O(3^n)\)，可以在更短时间内通过本题。

时间复杂度证明

设一个状态有 \(i\) 个元素，其子集一共有 \(2^i\) 个。有 \(i\) 个元素的状态一共有 \(C_n^i\) 种方案，总转移时间复杂度为：

\[ \sum_{i=1}^n{C_n^i\times 2^i}=3^n \]

此种方法比暴力转移要更优。

简单状压

Y556 最短路径

给定一个 \(n\) 个点 \(m\) 条边的有向图，有 \(k\) 个标记点，要求从规定的起点按任意顺序经过所有标记点到达规定的终点，问最短的距离是多少。

\(n\le 5\times 10^4,\ m\le 10^5,\ k \le 10\)。

我们可以在图上跑一遍 Dijkstra，同时通过状压记录经过了哪些标记点。时间复杂度 \(O((n\times 2^k + m\times 2 ^ k)\log (m\times 2^k))\)，不能通过本题，考虑优化。

通过仔细分析该算法，我们发现各个标记点之间的最短路被重复计算了 \(2^k\) 次，这是完全不必要的，考虑优化这些重复计算。

我们可以先求出这 \(k\) 个点之间两两的最短路，再在这 \(k\) 个点的简化图上跑状压 Dijkstra。这样，时间复杂度就被优化为 \(O\big(k(n+m)\log m+k^2\log k^2\big)\)，可以通过本题。

代码

合理利用 namespace 区分名称相似的函数可以避免很多不必要的麻烦。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
#include<iostream>
#include<cstring>
#include<queue>
#define int long long

using namespace std;
const int N = 5E4 + 10;
const int M = 1E5 + 10;
const int K = 15;

struct Edge {
    int v, w;
    int next;
} pool[2 * M];
int nn, head[N];

void addEdge(int u, int v, int w) {
    pool[++nn] = {v, w, head[u]};
    head[u] = nn;
}

int n, m, k, s, t;
int tag[N];
int b[K], dis[K][K];

namespace Dij1 {

    struct Node {
        int u, d;
        inline bool operator<(const Node &other) const {
            return d > other.d;
        }
    };

    int vis[N], d[N];
    priority_queue<Node> que;

    void Dijkstra(int x) {

        while(!que.empty()) que.pop();
        memset(d, 0x3f, sizeof(d));
        memset(vis, 0, sizeof(vis));
        d[b[x]] = 0;
        que.push({b[x], 0});
        while(!que.empty()) {
            int u = que.top().u;
            que.pop();
            if(vis[u]) continue;
            vis[u] = 1;
            for(int i = head[u]; i; i = pool[i].next) {
                int v = pool[i].v, w = pool[i].w;
                if(d[u] + w < d[v]) {
                    d[v] = d[u] + w;
                    que.push({v, d[v]});
                }
            }
        }
        for(int i = 1; i <= k + 2; i++) {
            dis[x][i] = d[b[i]];
        }

    }

}

namespace Dij2 {

    struct Node {
        int u, sta;
        int d;
        inline bool operator<(const Node &other) const {
            return d > other.d;
        }
    };

    int vis[K][1 << 12], d[K][1 << 12];
    priority_queue<Node> que;

    void Dijkstra() {

        if(s == t && k == 0) {
            cout << 0 << endl;
            exit(0);
        }

        memset(d, 0x3f, sizeof(d));
        d[k + 1][tag[s] ? (1 << tag[s] - 1) : 0] = 0;
        que.push({k + 1, tag[s] ? (1 << tag[s] - 1) : 0, 0});
        while(!que.empty()) {
            int u = que.top().u;
            int sta = que.top().sta;
            que.pop();
            if(vis[u][sta]) continue;
            vis[u][sta] = 1;
            if(u == k + 2 && (sta & ((1 << k) - 1)) == (1 << k) - 1) {
                cout << d[u][sta] << endl;
                exit(0);
            }
            for(int v = 1; v <= k + 2; v++) {
                if(u == v) continue;
                int tmp = sta | (1 << v - 1);
                if(d[u][sta] + dis[u][v] < d[v][tmp]) {
                    d[v][tmp] = d[u][sta] + dis[u][v];
                    que.push({v, tmp, d[v][tmp]});
                }
            }
        }

    }

}

signed main() {

    cin >> n >> m >> k >> s >> t;
    for(int i = 1; i <= m; i++) {
        int u, v, w;
        cin >> u >> v >> w;
        addEdge(u, v, w);
    }
    for(int i = 1; i <= k; i++) {
        int x;
        cin >> x;
        tag[x] = i;
        b[i] = x;
    }

    b[k + 1] = s;
    b[k + 2] = t;
    for(int i = 1; i <= k + 2; i++) {
        Dij1::Dijkstra(i);
    }

    Dij2::Dijkstra();

    cout << -1 << endl;

    return 0;
}

最小表示法

P2109 [NOI2007] 生成树计数

给定排成一行的 \(n\) 个节点，相邻节点之间的距离都为 \(1\)。每个点和所有与其距离不超过 \(k\) 的节点连边。问这个图的生成树的数量对 \(65521\) 取模的结果。

\(n\le 10^{15},k\le 5\)

\(n\le 10^{15}\)，但不易找到通项公式，因此考虑矩阵快速幂优化 dp。

考虑第 \(i\) 个节点，枚举其向 \([i-k,i-1]\) 中的哪些节点有连边。注意到 \(k\le 5\)，因此当前节点最多只能与前面的 \(5\) 个节点连边，dp 也只需要记录最近 \(5\) 个节点的状态。然而生成树要求图上不能出现环，因此我们需要保证第 \(i\) 个节点连接的 \([i-k,i-1]\) 里的点原本是不连通的。这就需要我们记录最后 \(k\) 个点的连通性。

要想记录所有可能出现的连通性状态，可以先暴力搜索出所有合法的连边方案，再把连通性相同的状态合并到一起，重新编号（这叫做最小表示法）。经过实测，\(k=5\) 时只有 \(52\) 种连通性状态。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
int id[55560], st[60], cnt[60], scnt;
int fa[K];

int find(int x) {
    if(fa[x] == x) return x;
    return fa[x] = find(fa[x]);
}

// 参数 s 表示边集 E 的二进制状压形式。
void check(int s) {
    for(int i = 1; i <= k; i++) fa[i] = i;
    for(int i = 1, c = 0; i <= k; i++) {
        for(int j = 1; j < i; j++, c++) {
            // 判环
            if(s & (1 << c)) {
                if(find(i) == find(j)) return;
                fa[max(find(i), find(j))] = min(find(i), find(j));
            }
        }
    }
    int sta = 0;
    for(int i = 1; i <= k; i++) {
        sta = sta * 10 + find(i);
    }
    if(!id[sta]) {
        id[sta] = ++scnt;
        st[scnt] = sta;
    }
    cnt[id[sta]]++;
}

int main() {

    // something...

    for(int i = 0; i < (1 << (k * (k - 1) / 2)); i++) {
        check(i);
    }

    // something...

}

这样，我们就得到了所有连通性状态的最小表示法和出现次数。接下来，我们只需要枚举第 \(i\) 个节点向 \([i-k,i-1]\) 的连边，判断是否有环，再找出新的连通性状态的最小表示法，就可以知道 \([i-k-1,i-2]\) 的每一个连通性状态 \(s_1\) 对 \([i-k,i-1]\) 的连通性状态 \(s_2\) 的贡献（转移系数）。

我们将转移系数写成一个矩阵，并将 \([i-k,i-1]\) 的每种连通性状态对应的方案数写在一个 \(1\times n\) 的向量里，将向量和矩阵相乘，即可得到后一个位置 \(i+1\) 的状态。

最后使用矩阵快速幂计算出 \(op^{n-k}\)，再以 \(cnt\) 为初始状态，得到最终状态。输出 \(k\) 个节点全部连通的状态对应的方案数即可。

代码 
#include<iostream>
#include<cstring>
#include<cmath>
#include<vector>
#define cint const int&
#define int long long
using namespace std;
const int K = 7;
const int MOD = 65521;
const int F[] = {-1, -1, 11, 111, 1111, 11111};

struct Matrix {
    int n;
    int a[60][60];
    inline Matrix(int _n) { n = _n; memset(a, 0, sizeof(a)); }
    inline int* operator [] (int index) { return a[index]; }
    inline const int* operator [] (int index) const { return a[index]; }
    inline Matrix operator * (const Matrix &b) const {
        Matrix res(n);
        for(int i = 1; i <= n; i++) {
            for(int j = 1; j <= n; j++) {
                for(int k = 1; k <= n; k++) {
                    res[i][j] = (res[i][j] + a[i][k] * b[k][j] % MOD) % MOD;
                }
            }
        }
        return res;
    }
};

inline vector<int> operator*(const vector<int> &a, const Matrix &b) {
    vector<int> res(b.n + 1);
    for(int i = 1; i <= b.n; i++) {
        for(int j = 1; j <= b.n; j++) {
            res[i] = (res[i] + a[j] * b[j][i]) % MOD;
        }
    }
    return res;
}

int n, k, lim;

int id[55560], st[60], cnt[60], scnt;
vector<int> beg, ed;

int p[K], fa[K];

int find(int x) {
    if(fa[x] == x) return x;
    return fa[x] = find(fa[x]);
}

void check(int s) {
    for(int i = 1; i <= k; i++) fa[i] = i;
    for(int i = 1, c = 0; i <= k; i++) {
        for(int j = 1; j < i; j++, c++) {
            if(s & (1 << c)) {
                if(find(i) == find(j)) return;
                fa[max(find(i), find(j))] = min(find(i), find(j));
            }
        }
    }
    int sta = 0;
    for(int i = 1; i <= k; i++) {
        sta = sta * 10 + find(i);
    }
    if(!id[sta]) {
        id[sta] = ++scnt;
        st[scnt] = sta;
    }
    cnt[id[sta]]++;
}

inline void add(Matrix &op, cint s1, cint s2) {
    op[id[s1]][id[s2]]++;
}

Matrix qpow(Matrix a, int k) {
    Matrix res(a.n);
    for(int i = 1; i <= a.n; i++) res[i][i] = 1;
    while(k) {
        if(k & 1) res = res * a;
        a = a * a;
        k >>= 1;
    }
    return res;
}

signed main() {

    cin >> k >> n;

    for(int i = 0; i < (1 << (k * (k - 1) / 2)); i++) {
        check(i);
    }

    beg.resize(scnt + 1);
    for(int i = 1; i <= scnt; i++) beg[i] = cnt[i];

    Matrix op(scnt);
    for(int s = 1; s <= scnt; s++) {
        int sta = st[s];
        for(int i = k; i >= 1; i--) p[i] = sta % 10, sta /= 10;
        for(int i = 0; i < (1 << k); i++) {
            for(int j = 1; j <= k; j++) fa[j] = p[j];
            fa[k + 1] = k + 1;
            bool flag = true;
            for(int j = 1; j <= k; j++) {
                if(i & (1 << j - 1)) {
                    if(find(k + 1) == find(j)) {
                        flag = false;
                        break;
                    }
                    fa[max(find(k + 1), find(j))] = min(find(k + 1), find(j));
                }
            }
            if(!flag) continue;
            int sta2 = 0, tag = 0;
            for(int i = 2; i <= k + 1; i++) if(find(i) == 1) { tag = i; break; }
            if(!tag) continue;
            for(int i = 2; i <= k + 1; i++) {
                sta2 = sta2 * 10 + (find(i) == 1 ? tag : find(i)) - 1;
            }
            add(op, st[s], sta2);
        }
    }

    op = qpow(op, n - k);
    ed = beg * op;

    cout << ed[id[F[k]]] << endl;

    return 0;
}

0420模拟赛 最小生成树

给定一张 \(n\) 个点 \(m\) 条边的 无向连通图，图中的第 \(i\) 条边有 \(a_i,b_i\)​ 两个权重。

对于图中的每条边，你可以选择 \(a_i,b_i\)​ 其中之一作为该条边的边权。

对于所有满足 \(0\le k\le m\) 的整数 \(k\)，你需要求出，若选择 恰好 \(k\) 个 \(a_i\)​ 作为对应边的边权，恰好 \(m−k\) 个 \(b_i\)​ 作为对应边的边权，该图的最小生成树的边权和 最大 是多少。

\(n\le 9,\ m\le 100\)

\(n\le 9\) 可以状压 dp。我们设 \(f[i][s][j]\) 表示考虑了前 \(i\) 条边，包含了 \(j\) 条 \(a\) 边，连通性为 \(s\)，最大的 mst 边权和是多少。

为了保证求得的是最小生成树，我们应该像 kk 算法一样，按边权从小到大考虑。对于一条边 \((u,v,a,b)\)，我们将它拆成两条边 \((u,v,a)\) 和 \((u,v,b)\)，然后将所有边按边权升序排序。

我们现在尝试向森林 \(s\) 中加入一条边 \((u,v,w)\)，分讨几种情况：

  • 森林中 \(u,v\) 未连通：
  • 是 \(a\) 边：
    • \(a\le b\)，可选可不选；
    • \(a>b\)，说明排在前面的 \(b\) 边没有选，由于 \(a,b\) 边必须要选一条，因此这条边必须选；
  • 是 \(b\) 边：
    • \(a\le b\)，说明排在前面的 \(a\) 边没有选，由于 \(a,b\) 边必须要选一条，因此这条边必须选；
    • \(a>b\)，可选可不选；
  • 森林中 \(u,v\) 已经连通：是否选择这条边不会对森林产生影响；

考虑这些情况如何进行转移。

  • 森林中 \(u,v\) 未连通：
  • 是 \(a\) 边：
    • \(a\le b\)，\(f[i-1][s][j]\to f[i][s'][j+1],\ f[i][s][j]\)
    • \(a>b\)，\(f[i-1][s][j]\to f[i][s'][j+1]\)
  • 是 \(b\) 边：
    • \(a>b\)，\(f[i-1][s][j]\to f[i][s'][j],\ f[i][s][j]\)
    • \(a\le b\)，\(f[i-1][s][j]\to f[i][s'][j]\)

然而对于森林中 \(u,v\) 已经连通的情况，我们不知道连接 \(u,v\) 两个连通块的是 \((u,v)\) 对应的另一条边，还是一条不是 \((u,v)\) 的边。因此我们不知道对应的另一条边有没有被选择，也就不知道这一条边是否应该选择，不易处理 \(j\) 的转移。

考虑将 \(j\) 的贡献放在边权较小，先被考虑的一条边上面统计：

  • 森林中 \(u,v\) 未连通：
  • 是 \(b\) 边且 \(a>b\) 的情况：\(f[i-1][s][j]\to f[i][s'][j],\ f[i][s][j+1]\)
  • 森林中 \(u,v\) 已连通：
  • 是 \(a\) 边且 \(a\le b\)，\(f[i-1][s][j]\to f[i][s][j],\ f[i][s][j+1]\)
  • 是 \(b\) 边且 \(a>b\)，\(f[i-1][s][j]\to f[i][s][j+1],\ f[i][s][j]\)
  • 其余情况，\(f[i-1][s][j]\to f[i][s][j]\)，因为贡献已经被另一条边计算过了；

我们搜索出所有连通性状态 \(s\)，在 \(n=9\) 时有 \(Bell(9)=21147\) 种状态。时间复杂度为 \(O\big(m^2Bell(9)\big)\)。

由于 \(n\le 9\)，我们不能使用 \(2^{\binom n2}\) 的纯暴力枚举。考虑枚举当前节点向之前的哪些节点连边，再通过哈希使得每层的每一种连通性只会被处理一遍，可以将时间复杂度降低到 \(O\big(n2^nBell(n)\big)\)，可以通过本题，并且实际跑不满。

代码 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#include<iostream>
#include<algorithm>
#include<unordered_set>
#include<unordered_map>
#include<cstring>
using namespace std;
const int N = 10;
const int M = 110;
const int INF = 0x3f3f3f3f;

#define gmax(x, y) x = max(x, y)

struct Edge {
    int u, v, w, t;
    inline bool operator<(const Edge &b) const {
        if(w != b.w) return w < b.w;
        return t > b.t;
    }
} edg[2 * M];
int ne;

int n, m;
int conn[22000][N], nn, ed;

int fa[N][N];
int find(int id, int x) { if(fa[id][x] == x) return x; return fa[id][x] = find(id, fa[id][x]); }
void merge(int id, int x, int y) { x = find(id, x); y = find(id, y); if(x < y) swap(x, y); fa[id][x] = y; }
unordered_set<int> his[N];
unordered_map<int, int> stat;

int hsh(int id) {
    int res = 0;
    for(int i = 0; i < n; i++) find(id, i);
    for(int i = 0; i < n; i++) res = res * n + fa[id][i];
    return res;
}

void dfs(int x) {
    int hs = hsh(x);
    if(his[x].count(hs)) return;
    his[x].insert(hs);
    if(x == n) {
        ++nn;
        for(int i = 0; i < n; i++) conn[nn][i] = fa[x][i];
        stat[hs] = nn;
        bool flag = true;
        for(int i = 0; i < n; i++) flag &= (fa[x][i] == 0);
        if(flag) ed = nn;
        return;
    }
    for(int i = 0; i < (1 << x); i++) {
        for(int j = 0; j < n; j++) fa[x + 1][j] = fa[x][j];
        for(int j = 0; j < x; j++) {
            if(i & (1 << j)) {
                merge(x + 1, x, j);
            }
        }
        dfs(x + 1);
    }
}

int f[22000][M], g[22000][M];

// #define FIO

int main() {

    #ifdef FIO
    freopen("mst.in", "r", stdin);
    freopen("mst.out", "w", stdout);
    #endif

    cin >> n >> m;
    for(int i = 1; i <= m; i++) {
        int u, v, a, b;
        cin >> u >> v >> a >> b;
        --u, --v;
        if(a <= b) {
            edg[++ne] = {u, v, a, 1};
            edg[++ne] = {u, v, b, 0};
        } else {
            edg[++ne] = {u, v, b, -1};
            edg[++ne] = {u, v, a, 0};
        }
    }

    sort(edg + 1, edg + 1 + ne);

    for(int i = 0; i < n; i++) fa[1][i] = i;
    dfs(1);

    for(int s = 0; s <= nn; s++) {
        for(int j = 0; j <= ne; j++) {
            g[s][j] = -INF;
        }
    }
    g[1][0] = 0;

    for(int i = 0; i < ne; i++) {
        int u = edg[i + 1].u, v = edg[i + 1].v, w = edg[i + 1].w, t = edg[i + 1].t;
        for(int s = 1; s <= nn; s++) {
            for(int j = 0; j <= m; j++) {
                f[s][j] = g[s][j];
                g[s][j] = -INF;
            }
        }
        for(int s = 1; s <= nn; s++) {
            for(int k = 0; k < n; k++) fa[0][k] = conn[s][k];
            merge(0, u, v);
            int s1 = stat[hsh(0)], flag = (s1 != s);
            for(int j = 0; j <= m; j++) {
                if(t == 1) {
                    gmax(g[s1][j + 1], f[s][j] + w * flag);
                    gmax(g[s][j], f[s][j]);
                } else if(t == -1) {
                    gmax(g[s1][j], f[s][j] + w * flag);
                    gmax(g[s][j + 1], f[s][j]);
                } else {
                    gmax(g[s1][j], f[s][j] + w * flag);
                }
            }
        }
    }

    for(int i = 0; i <= m; i++) {
        cout << g[ed][i] << '\n';
    }

    return 0;
}

多行状态

P2704 [NOI2001] 炮兵阵地

YbtOj 链接

题面请参考洛谷链接。

我们发现第 \(i\) 行的状态受到第 \(i-1\) 行和 \(i-2\) 行的影响。状态应当同时包含 \(i\) 行和 \(i-1\) 行的信息。

我们设状态 \((i,j_1,j_2)\) 表示第 \(i\) 行的状态为 \(j_1\)，第 \(i-1\) 行的状态为 \(j_2\)，最多能摆放多少个大炮。

\[ f_{i,j_1,j_2}=\max_{j_3\cap j_1=\emptyset\text{ and }j_3\cap j_2=\emptyset}\{f_{i-1,j_2,j_3}+cnt(j_1)\} \]

再通过预处理每行的 valid 数组，找到最小表示法，能极大优化时间复杂度（因为每行的合法状态不超过 \(70\) 个）。

多维状压

有时，单个点的状态数可能大于 \(2\)，此时无法用二进制简单的表示出此状态。为此，对于一个状态 \(i\)，其第 \(j\) 个位置的状态可以表示为 \(\lfloor\frac{i}{k^j}\rfloor \bmod k\)。这样，对于 \(k\) 维的信息，我们也能使用状压处理了。只不过此种方法无法像二进制一样可以进行位运算。

P8504 「CGOI-2」No will to break

一场战斗由 \(n\) 个时刻组成，第 \(i\) 个时刻有 \(\frac{x_i}{x_i+y_i}\) 的概率是安全的。

在安全的时刻，你可以进行“聚集”。要求每连续的 \(a\) 个时刻都至少要有 \(b\) 个时刻进行聚集，在此前提下希望进行聚集的时刻数量尽量少；若不能满足此前提则认为进行聚集的时刻数量为 \(0\)。求进行聚集的时刻数量的期望，答案对 \(998244353\) 取模。

\(n\le 1.5\times 10^4,~1\le b<a\le\min(n,9)\)

在已知哪些点是安全的情况下，显然有贪心策略：能不放就不放，尽可能往后放。而又注意到我们对一个前缀只关心其后 \(a-1\) 位，因此考虑状压这些位置。发现一个位置只有三种状态：安全不聚集、安全聚集、不安全，因此考虑用三进制表示状态。

考虑在一个 \(a-1\) 位状态后面添加一位 \(i\)，如何转移。如果 \([i-a+1,i]\) 中聚集的数量不超过 \(b\) 个，我们就从右往左依次考虑每一个安全不聚集的位置，进行一次聚集（\(i\) 初始是不聚集的），直到聚集的次数达到 \(b\) 我们就得到了转移后的状态，或者不存在安全不聚集的位置但次数仍然不足 \(b\)，此时转移不合法。

考虑如何转移期望，以及如何处理概率。记 \(f_{i,S}\) 表示（钦定）前缀 \(i\) 的后 \(a-1\) 位状态为 \(S\) 的前提下，聚集次数的期望（条件期望）（不合法算 \(0\)）。注意到一次转移可能会给聚集次数 \(+1\) 或者不加，然而 \(+1\) 只对原本合法的状态有贡献（原本不合法就还是 \(0\)），因此加的 \(1\) 需要乘以原本合法的概率。设 \(g_{i,S}\) 表示后 \(a-1\) 位被钦定的前提下，合法的概率，也容易转移。同时往后挪动一位会导致 \(i-a+1\) 从钦定变为不钦定，因此再乘以这一位安全/不安全的概率。转移如下：

\[ f_{i,\delta(j,c)}\gets (f_{i-1,j}+g_{i-1,j}w_{j,c})p_{i-a+1,high(j)}\\ g_{i,\delta(j,c)}\gets g_{i-1,j}p_{i-a+1,high(j)} \]

\(high(j)\) 表示状态 \(j\) 的最高位（\(i-a+1\)）是否安全。从 \(a\) 开始转移，先把前 \(a\) 个一次性聚集 \(b\) 次。统计答案时就枚举 \(f_n\) 的所有状态，把钦定的后 \(a-1\) 位都乘以对应的概率转化为不钦定贡献到答案。这样转移可以做到 \(O(3^an)\)，不能通过。

发现除去两边的边界，中间有效的状态要么包含 \(b-1\) 次聚集要么包含 \(b\) 次聚集，否则要么不合法要么不符合贪心。因此中间转移时只考虑包含 \(b-1\) 或 \(b\) 次聚集的状态即可，差不多是一个组合数，数一下最大约 \(3584\)。\(high\) 要预处理，要不然开销很大。

代码 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#include<iostream>
#define ll long long
using namespace std;
const int N = 1.5e4 + 10;
const int K = 9;
const int MOD = 998244353;

inline void add(int &a, int b) { a = (a + b) % MOD; }

inline int qpow(int a, int b) {
    int res = 1;
    while(b) {
        if(b & 1) res = (ll)res * a % MOD;
        a = (ll)a * a % MOD;
        b >>= 1;
    }
    return res;
}

const int pw3[10] = {1, 3, 9, 27, 81, 243, 729, 2187, 6561, 19683};

int n, a, b, m;
int x[N][2];
int trans[6570][2], w[6570][2];

int cnt[6750];
int f[2][6570], g[2][6570];

int sta[6750], nn;
int high[6750];

int main() {

    cin >> n >> a >> b;
    for(int i = 1; i <= n; i++) cin >> x[i][0] >> x[i][1];

    m = pw3[a - 1];
    for(int i = 0; i < m; i++) {
        int s[K + 2] = {}, tc = 0;
        for(int j = 0; j < a - 1; j++) s[j + 1] = i / pw3[j] % 3, cnt[i] += (s[j + 1] & 1);
        if(cnt[i] >= b) { trans[i][0] = i * 3 % m, trans[i][1] = i * 3 % m + 2; continue; }
        tc = b - cnt[i];
        for(int c = 0; c <= 1; c++) {
            int v = i * 3 % m + (c << 1), p = -1;
            s[0] = c << 1;
            for(int k = 1; k <= tc; k++) {
                int j = p + 1; p = -1;
                for(; j < a; j++) {
                    if(s[j] == 0) { p = j; break; }
                }
                if(p == -1) break;
                v += pw3[p];
            }
            if(p == -1) trans[i][c] = -1;
            else trans[i][c] = v % m, w[i][c] = tc;
        }
    }
    for(int i = 0; i < m; i++) {
        if(cnt[i] == b || cnt[i] == b - 1) sta[++nn] = i;
    }

    for(int i = 1; i <= n; i++) {
        int d = qpow((x[i][0] + x[i][1]) % MOD, MOD - 2);
        x[i][0] = (ll)x[i][0] * d % MOD;
        x[i][1] = (ll)x[i][1] * d % MOD;
    }

    for(int j = 0; j < m; j++) {
        bool flag = true;
        for(int i = 0; i < a; i++) {
            if(j / pw3[i] % 3 == 1) { flag = false; break; }
        }
        g[(a & 1) ^ 1][j] = flag;
    }
    --a;

    for(int j = 0; j < m; j++) high[j] = ((j << 1) + j) / m >> 1;

    for(int j = 0; j < m; j++) f[(a + 1) & 1][j] = g[(a + 1) & 1][j] = 0;
    for(int i = a + 1; i <= n; i++) {
        int t = i & 1;
        for(int j = 1; j <= nn; j++) f[t][sta[j]] = g[t][sta[j]] = 0;
        if(i == a + 1)
            for(int j = 0; j < m; j++) {
                int c = 0;
                if(~trans[j][c]) {
                    add(f[t][trans[j][c]], (f[t ^ 1][j] + (ll)w[j][c] * g[t ^ 1][j] % MOD) * x[i - a][high[j]] % MOD);
                    add(g[t][trans[j][c]], (ll)g[t ^ 1][j] * (x[i - a][high[j]]) % MOD);
                }
                c = 1;
                if(~trans[j][c]) {
                    add(f[t][trans[j][c]], (f[t ^ 1][j] + (ll)w[j][c] * g[t ^ 1][j] % MOD) * x[i - a][high[j]] % MOD);
                    add(g[t][trans[j][c]], (ll)g[t ^ 1][j] * (x[i - a][high[j]]) % MOD);
                }
            }
        else
            for(int jj = 1; jj <= nn; jj++) {
                int j = sta[jj];
                int c = 0;
                if(~trans[j][c]) {
                    add(f[t][trans[j][c]], (f[t ^ 1][j] + (ll)w[j][c] * g[t ^ 1][j] % MOD) * x[i - a][high[j]] % MOD);
                    add(g[t][trans[j][c]], (ll)g[t ^ 1][j] * (x[i - a][high[j]]) % MOD);
                }
                c = 1;
                if(~trans[j][c]) {
                    add(f[t][trans[j][c]], (f[t ^ 1][j] + (ll)w[j][c] * g[t ^ 1][j] % MOD) * x[i - a][high[j]] % MOD);
                    add(g[t][trans[j][c]], (ll)g[t ^ 1][j] * (x[i - a][high[j]]) % MOD);
                }
            }
    }

    int ans = 0;
    for(int i = 0; i < m; i++) {
        int res = f[n & 1][i];
        for(int j = 1; j <= a; j++) res = (ll)res * x[n - j + 1][(i / pw3[j - 1]) % 3 >> 1] % MOD;
        ans = (ans + res) % MOD;
    }
    cout << ans << '\n';

    return 0;
}