NOI2015 寿司晚宴

在晚宴上，主办方为大家提供了 n−1 种不同的寿司，编号 1,2,3,\ldots,n-1，其中第种寿司的美味度为 i+1。（即寿司的美味度为从 2 到 n）

现在小 G 和小 W 希望每人选一些寿司种类来品尝，他们规定一种品尝方案为不和谐的当且仅当：小 G 品尝的寿司种类中存在一种美味度为 x 的寿司，小 W 品尝的寿司中存在一种美味度为 y 的寿司，而 x 与 y 不互质。

现在小 G 和小 W 希望统计一共有多少种和谐的品尝寿司的方案（对给定的正整数 p 取模）。注意一个人可以不吃任何寿司。

n \leq 300,p \leq 10^{10}

解题思路：

状压 DP。

对于 30 分的部分分，可以考虑设计一个暴力 DP。这里使用一个小 trick，我们只需要记录每个数的质因子即可，而在 30 以内的质数一共只有 10 个，故我们不妨设计 f_{i,s_1,s_2} 表示选到第 i 个数，小 G 的寿司质因子集合为 s_1，小 W 的寿司质因子集合为 s_2 ，这个可以压成一个二进制数来存。

那么转移的时候我们就考虑 i 放到 s_1 与 s_2 中的方案，不难写出转移方程（其中 k 是 i 的质因数二进制表示的数）：

f_{i,s_1 \text{ or } k ,s_2} = f_{i,s_1 \text{ or } k ,s_2} + f_{i-1,s_1,s_2} \text{ (if s1 and k = 0)} \\ f_{i,s_1 ,s_2 \text{ or } k} = f_{i,s_1,s_2 \text{ or } k} + f_{i-1,s_1,s_2} \text{ (if s2 and k = 0)} \\

这个转移是比较简单的，时间复杂度 \mathcal{O}(n \times {2^{10}}^2)，注意这个第一维是完全可以滚掉的，滚掉的方法可以参考数组，这边也是记录一下这个转移的方法。

接下来考虑满分做法，我们发现瓶颈在于我们直接把所有质因子压进了状态，这样在 n \leq 500 时最多会有 95 个因子，无法接受。

这就要用到本题启发的第二个 trick：对于具有特殊性的一些因子单独处理。具体来说，我们发现所有 \geq \sqrt{n} 的质因子最多只会出现一次，证明很显然，如果这个因子 p \leq \sqrt{n}，那么有 p^2 \leq n ，那么这样我们状态里只要压缩 \leq \sqrt{n} = 22 的所有质数，一共只有八个。

所以对于 \geq \sqrt{n} 的因子我们就要单独考虑了，对于一个质因子 p ，我们首先考虑都能表示成 kp(k\in \N^*,k < p) 的形式，然后我们考虑对于所有的 kp 来说，我们需要单独计算，也就是说我们把这个数列按照 p 来分段。

  • 如果 p 相同，那么我们单独考虑把这个数分给 s_1,s_2 的方案数，但是这个结果不能并入 f 数组中，因为对于后面的 kp 会造成影响。所以我们不妨新建两个数组 g,h 来表示把该数分入 s_1,s_2 的方案数。

    对于 g 来说（对于 h 是类似的），我们考虑一共两种情况：

    1. 在已经有 p 的基础上加入这个数，那么这部分的方案就是 g_{i-1,s_1,s_2}
    2. 在没有 p 的基础上加入这个数，这部分的方案就是 f_{i,s_1,s_2}

    个人认为这种思路是比题解里大部分的容斥更加自然的。

  • 如果 p 不同，我们直接继承上一段的状态，也就是 f 正常更新，把 g,h 赋值为 f 即可。

这样的时间复杂度是很优秀的！达到了 \mathcal{O}(n \times {2^8}^2)

代码：

const int N = 510;
const int px[8] = {2,3,5,7,11,13,17,19};

int f[2][(1<<8)][(1<<8)],g[2][(1<<8)][(1<<8)],w[2][(1<<8)][(1<<8)],mod,n,tot,prim[N],to[500],has[N],lis[N],ans;
bool vis[31];
vector <int> fac[N];

bool cmp(int a,int b) {
    if (prim[a] != prim[b]) return prim[a] < prim[b];
    return a / prim[a] < b / prim[b];
}

signed main() {
    read(n,mod);
    for (int i = 2;i <= n;++i) {
        int x = i;
        for (int j = 0;j < 8;++j) {
            if (x % px[j] == 0) {
                has[i] ^= (1 << j);
                while (x % px[j] == 0) x /= px[j];
            }
        }
        prim[i] = x;
        lis[i] = i;
    }
    sort(lis + 2,lis + 1 + n,cmp);
    f[1][0][0] = 1;
    int cur = 1,lst = 0;
    for (int i = 2;i <= n;++i) {
        int x = lis[i];
        //printf("%d %d\n",x,prim[x]);
        cur ^= 1,lst ^= 1;
        for (int s1 = 0;s1 < (1 << 8);++s1) 
            for (int s2 = 0;s2 < (1 << 8);++s2) 
                f[cur][s1][s2] = f[lst][s1][s2],w[cur][s1][s2] = w[lst][s1][s2],g[cur][s1][s2] = g[lst][s1][s2];
        for (int s1 = 0;s1 < (1 << 8);++s1) {
            for (int s2 = 0;s2 < (1 << 8);++s2) {
                if (s1 & s2) continue;
                if (prim[x] == 1) {//nobigprime
                    if ((s2 & has[x]) == 0) f[cur][s1 | has[x]][s2] = (f[cur][s1 | has[x]][s2] + f[lst][s1][s2]) % mod;
                    if ((s1 & has[x]) == 0) f[cur][s1][s2 | has[x]] = (f[cur][s1][s2 | has[x]] + f[lst][s1][s2]) % mod;
                } else if (prim[x] != prim[lis[i - 1]]) {
                    f[cur][s1][s2] = (f[lst][s1][s2] + w[lst][s1][s2] + g[lst][s1][s2]) % mod;
                    g[cur][s1][s2] = w[cur][s1][s2] = f[cur][s1][s2];
                } else {
                    if ((s2 & has[x]) == 0) g[cur][s1 | has[x]][s2] = ((g[cur][s1 | has[x]][s2] + g[lst][s1][s2]) % mod+ f[cur][s1][s2]) % mod;
                    if ((s1 & has[x]) == 0) w[cur][s1][s2 | has[x]] = ((w[cur][s1][s2 | has[x]] + w[lst][s1][s2]) % mod + f[cur][s1][s2]) % mod;
                }
            }
        }
    }
    for (int s1 = 0;s1 < (1 << 8);++s1) {
        for (int s2 = 0;s2 < (1 << 8);++s2) {
            if ((s1 & s2) == 0) ans = (ans + (f[cur][s1][s2]+ g[cur][s1][s2]) % mod + w[cur][s1][s2]) % mod;
        }
    }
    printf("%lld\n",ans);
    return 0;
}