[HAOI2017] 方案数 题解

教练布置下来当容斥单元小测的，磕了两个小时终于磕出来了，纪念一下。

其实题目应该说明一下每次操作后坐标不能与原坐标相同。

---

首先不考虑障碍物的限制，如何求出从 $(0, 0, 0)$ 走到 $(x, y, z)$ 的方案数？

我们很容易想出一个非常暴力的dp：$dp_{i, j, k}$ 表示坐标为 $(i, j, k)$ 时从 $(0, 0, 0)$ 走到当前位置的方案数，直接枚举当前位置的上一个可能的位置来转移。

但是这样的时间复杂度太大了。

观察题目条件，假设我们某一步是把 $x$ 更改为 $x'$，因为 $x \land x' = x$，所以二进制下 $x'$ 的 $1$ 的数量一定是大于 $x$ 的，这启发我们用二进制下 $1$ 的个数来作为dp的状态。

重新定义 $dp_{i, j, k}$ 表示满足 $\operatorname{popcount}(x) = i, \operatorname{popcount}(y) = j, \operatorname{popcount}(z) = k$ 的坐标 $(x, y, z)$ 的方案数，其中 $\operatorname{popcount}(x)$ 表示二进制下 $x$ 含有的数字 $1$ 的个数。

因为每次转移中，有且仅有某一维的状态改变，并且是增大，所以只需要考虑这一个增加量有多少种可能即可。

假设这一增加量为 $l$，从 $dp_{i - l, j, k}$ 转移到 $dp_{i, j, k}$，那么这一增加量 $l$ 相当于是在 $i$ 个二进制位中选取 $l$ 个二进制位，可能的数量为 $\binom{i}{l}$。

如果要再说细一点，就是因为每一个二进制位相互之间有差异，但是任意两个二进制位被选取的先后对答案没有影响，所以是用组合数计算而不是排列数或者直接认为是 $1$。

状态转移方程就可以写出来了：

$$ dp_{i, j, k} = \sum\limits_{l = 1}^{i}\binom{i}{l}dp_{i - l, j, k} + \sum\limits_{l = 1}^{j}\binom{j}{l}dp_{i, j - l, k} + \sum\limits_{l = 1}^{k}\binom{k}{l}dp_{i, j, k - l} $$

边界值：$dp_{0, 0, 0} = 1$。

为什么这样子dp过后坐标为 $(x, y, z)$ 的方案数就是 $dp_{\operatorname{popcount}(x), \operatorname{popcount}(y), \operatorname{popcount}(z)}$ 了呢，不会多统计到其它的路径吗？

不会，因为我们在dp的过程中是对当前状态的二进制下有 $1$ 的每一位考虑的，换言之，我们是每次若干位地填充这个状态的，并且不会多填。

接下来考虑有障碍的情况该怎么做。

（如果做过“两双手”这道题的同学应该能直接看出来出来这是dp+容斥。）

首先要把那些 $x \land n \not = x$ 或 $y \land m \not = y$ 或 $z \land r \not = z$ 的坐标排除掉，因为从这个坐标不可能走到终点。

根据上面的结论：每次移动二进制下 $1$ 的个数只增不减，我们可以按照 $x, y , z$ 的二进制下 $1$ 的个数为优先级将这些坐标排序，这样影响某个位置的方案的只可能是数组中在它前面的元素。

我们发现如果要计算恰好所有的障碍点都不经过的方案会很麻烦，考虑将恰好这一要求转化一下。

我们可以统计最后一步走完后停留在某个障碍点并且之前没有经过其它障碍点的方案数，相当于把恰好不经过所有的障碍点转化为了不经过前面若干个障碍点，这也是容斥的精髓。

于是令 $f_{i}$ 表示停留在第 $i$ 个障碍点并且不经过第 $1 \sim i - 1$ 个障碍点的方案数，可以得到：

（下面将 $\operatorname{popcount}(x)$ 简写为 $\operatorname{cnt}(x)$）

$$ f_{i} = dp_{\operatorname{cnt}(x_{i}), \operatorname{cnt}(y_{i}), \operatorname{cnt}(z_{i})} - \sum\limits_{j = 1}^{i - 1}f_{j}dp_{\operatorname{cnt}(x_{i} \oplus x_{j}), \operatorname{cnt}(y_{i} \oplus y_{j}), \operatorname{cnt}(z_{i} \oplus z_{j})} $$

上面的 $\oplus$ 表示按位异或。

注意这里枚举的第 $j$ 个障碍点必须满足其能够到达第 $i$ 个障碍点，也就是 $x_{j} \land x_{i} = x_{j}, y_{j} \land y_{i} = y_{j}, z_{j} \land z_{i} = z_{j}$，这一段式子太长了就没写在上面。

解释一下上面的式子：不考虑前面是否经过别的障碍点的方案数是 $dp_{\operatorname{cnt}(x_{i}), \operatorname{cnt}(y_{i}), \operatorname{cnt}(z_{i})}$，然后需要减去经过了障碍点的方案数，于是枚举经过的第一个障碍点 $j$，而第 $j$ 个障碍点到第 $i$ 个障碍点的路径方案数为 $dp_{\operatorname{cnt}(x_{i} \oplus x_{j}), \operatorname{cnt}(y_{i} \oplus y_{j}), \operatorname{cnt}(z_{i} \oplus z_{j})}$，因为二进制下已经为 $1$ 的点不用再考虑了，所以要考虑的是坐标按位异或后的值，也就是二进制下相异的位数。

为什么可以只用考虑前面的障碍点而不考虑后面的？注意到我们是按照 $\operatorname{popcount}(x), \operatorname{popcount}(y), \operatorname{popcount}(z)$ 为优先级排序了的，在一个障碍点后面的点不可能影响到这个障碍点。

如何计算最后的答案？其实只需要把 $(n, m, r)$ 也当作一个障碍点丢进数组里就好了。

要注意 __builtin_popcount 传参传的是 unsigned int，要用 __builtin_popcountll！我被这个坑了。

附一份丑陋的代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 998244353;
ll n, m, r, o, x, y, z, cnt, fct[10005], inv[10005], dp[10005], f[64][64][64];
ll popcount(ll x) {
    ll ret = 0;
    while(x) ++ret, x ^= (x & (-x));
    return ret;
}
struct cube {
    ll x, y, z, wx, wy, wz, w;
    cube(ll X = 0, ll Y = 0, ll Z = 0): x(X), y(Y), z(Z) {
        wx = popcount(X);
        wy = popcount(Y);
        wz = popcount(Z);
        w = wx + wy + wz;
    }
    cube operator = (const cube& _) {
        x = _.x, y = _.y, z = _.z;
        wx = _.wx, wy = _.wy, wz = _.wz;
        w = _.w;
        return *this;
    }
    bool operator < (const cube& _) const {
        if(wx != _.wx) return wx < _.wx;
        if(wy != _.wy) return wy < _.wy;
        if(wz != _.wz) return wz < _.wz;
        return false;
    }
} c[10005];
ll ksm(ll x, ll y) {
    ll ret = 1;
    while(y) {
        if(y & 1) ret = ret * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ret;
}
ll C(ll n, ll m) {
    return fct[n] * inv[m] % mod * inv[n - m] % mod;
}
int main() {
    ios::sync_with_stdio(0);
    cin.tie(0), cout.tie(0);
    fct[0] = 1;
    for(int i = 1; i <= 10000; ++i) fct[i] = fct[i - 1] * i % mod;
    inv[10000] = ksm(fct[10000], mod - 2);
    for(int i = 10000; i >= 1; --i) inv[i - 1] = inv[i] * i % mod;
    cin >> n >> m >> r >> o;
    while(o--) {
        cin >> x >> y >> z;
        if((n & x) == x && (m & y) == y && (r & z) == z) {
            c[++cnt] = cube(x, y, z);
        }
    }
    c[++cnt] = cube(n, m, r);
    stable_sort(c + 1, c + 1 + cnt);
    f[0][0][0] = 1;
    for(int i = 0; i <= 63; ++i) {
        for(int j = 0; j <= 63; ++j) {
            for(int k = 0; k <= 63; ++k) {
                for(int l = 1; l <= i || l <= j || l <= k; ++l) {
                    if(i >= l) f[i][j][k] = (f[i][j][k] + C(i, l) * f[i - l][j][k] % mod) % mod;
                    if(j >= l) f[i][j][k] = (f[i][j][k] + C(j, l) * f[i][j - l][k] % mod) % mod;
                    if(k >= l) f[i][j][k] = (f[i][j][k] + C(k, l) * f[i][j][k - l] % mod) % mod;
                }
            }
        }
    }
    for(int i = 1; i <= cnt; ++i) {
        // cerr << i << ":\n";
        for(int k = 1; k < i; ++k) {
            if((c[i].x & c[k].x) == c[k].x && (c[i].y & c[k].y) == c[k].y && (c[i].z & c[k].z) == c[k].z) {
                // cerr << k << " ";
                dp[i] = (dp[i] - dp[k] * f[c[i].wx - c[k].wx][c[i].wy - c[k].wy][c[i].wz - c[k].wz] % mod + mod) % mod;
            }
        }
        dp[i] = (dp[i] + f[c[i].wx][c[i].wy][c[i].wz] % mod) % mod;
        // cerr << "\n> " << dp[i] << '\n';
    }
    cout << dp[cnt];
    return 0;
}