洛谷 P2150 - 寿司晚宴
我记得这是我 NOI 的某天晚上去自习室写的题，然后题解看不懂然后没电了，回寝室之后 AC 了，AC 了然后还是不会（
还有，洛谷咋这么喜欢关题解入口啊？？？
洛谷题目页面传送门
给定 \(n\)，求有多少个集合有序二元组 \((S,T)\) 满足 \(S,T\subseteq\{2,3,\cdots,n\}\) 且 \(\forall x\in S,\forall y\in T,x\perp y\)。答案对 \(p\) 取模。
\(n\in[2,500],p\in\left[1,10^{10}\right]\)。
显然，不能打表
显然，条件可以转化为：\(S,T\) 分别的质因数集合没有交集。于是很容易想到状压 DP。
容易想到 \(2\) 个比较暴力的状压 DP：
- \(dp_{i,j}\) 表示考虑到数 \(i\)，\(S\) 的质因数集合为 \(j\) 的方案数。很好转移，刷表即可；比较难的地方在于最后统计答案，注意到 \(S,T\) 的情况应该是一样的，所以答案就是 \(\sum\limits_{A\cap B=\varnothing}dp_{n,A}dp_{n,B}\)。时间复杂度 \(\mathrm O\!\left(n2^{\pi(n)}\right)\)。
- 不像上面那样绕弯子了，直白一点：\(dp_{i,j,k}\) 表示考虑到数 \(i\)，\(S,T\) 的质因数集合分别为 \(j,k\) 的方案数。转移依然刷表；统计答案就直接统计了。时间复杂度 \(\mathrm O\!\left(n4^{\pi(n)}\right)\)。
写个程序随便算一下发现 \(n=500\) 时 \(\pi(n)=95\)，无论哪种都是跑不过的。
注意到一个性质：一个数的所有质因数之积要小于等于原数。那么显然，\(\geq\sqrt x\) 的 \(x\) 的质因数只有一个。于是我们可以利用根号分治的基本思想，将每个数的质因数分为两类：小质因数和大质因数。显然大质因数的数量在 \(0\) 到 \(1\) 之间。把有大质因数的数，大质因数相同的分为一组；没有大质因数的数，每个数单独分为一组。这样显然，大质因数的限制仅在于每组之内，即每组最多只能有 \(1\) 个集合选数；然后带着这个限制，我们可以抛开大质因数只考虑小质因数了，考虑每组分别算，之后合并。小质因数最多只有 \(\pi(\sqrt n)=8\) 个，感觉很行。
考虑在小质因数范围内套用上面的暴力状压 DP。第 \(1\) 种是萎掉了，因为最后统计答案的时候还是要考虑大质因数，而我们没有把它们包含在状态里（所以你别看它复杂度小，其实局限性大。很多其他题也是这样的，标算往往是从比较 naive 的暴力优化过来的）；于是只能用第 \(2\) 种。加上这里大质因数带来的特殊限制条件，原来转移的时候有 \(3\) 种贡献方式：不分、分 \(S\)、分 \(T\)。而现在不行了，可以重新设 \(2\) 个 DP 数组 \(dp1,dp2\)，然后每个的转移只剩 \(2\) 种贡献方式。
接下来考虑如何合并两组。如果直接用两组 DP 数组最终状态直接进行合并的话，是四次方的复杂度，没有前途。不妨换一个思路，第二组从初始化 DP 数组的时候就把第一组的最终状态继承过来，最后随便容斥一下即可。这样对于每个数，都要有一个 \(\mathrm O\!\left(4^{\pi(\sqrt n)}\right)\) 的 DP 数组，毛估估一下理论可过，但是 \(\bmod\) 运算写的丑的话会被卡常，其实可以优化。空间可以滚动数组滚成总 \(\mathrm O\!\left(4^{\pi(\sqrt n)}\right)\) 的；时间上的话，注意到若 DP 数组两个维度有交集的话，一定是不合法的，直接不看，可以优化到 \(\mathrm O\!\left(n3^{\pi(\sqrt n)}\right)\)，但我比较懒，初始化的时候直接 memset 理论上复杂度没优化，但是大常数转移的时候优化了，可以通过。
代码：
#include<bits/stdc++.h>
using namespace std;
#define int long long
#define pb push_back
const int N=500;
int n,mod;
int dp[2][1<<8][1<<8];
int dp1[2][1<<8][1<<8],dp2[2][1<<8][1<<8];
vector<int> hav[N+1];
int id[N+1];
int mkmsk(vector<int> v,bool nobk=false){
int res=0;
for(int i=0;i+nobk<v.size();i++)res|=1<<id[v[i]];
return res;
}
signed main(){
id[2]=0;id[3]=1;id[5]=2;id[7]=3;id[11]=4;id[13]=5;id[17]=6;id[19]=7;
cin>>n>>mod;
for(int i=2;i<=n;i++){
vector<int> v;
int cpy=i;
for(int j=2;j*j<=cpy;j++)if(cpy%j==0){
v.pb(j);
while(cpy%j==0)cpy/=j;
}
if(cpy>1)v.pb(cpy);
if(v.size()&&v.back()>22)hav[v.back()].pb(mkmsk(v,true));
else hav[i].pb(mkmsk(v));
}
int now=0;
dp[0][0][0]=1;
for(int i=1;i<=n;i++)if(hav[i].size()){
now++;
for(int j=0;j<1<<8;j++)for(int k=0;k<1<<8;k++)dp1[0][j][k]=dp2[0][j][k]=dp[now-1&1][j][k];
for(int j=0;j<hav[i].size();j++){
memset(dp1[j+1&1],0,sizeof(dp1[j+1&1])),memset(dp2[j+1&1],0,sizeof(dp2[j+1&1]));
for(int k=0;k<1<<8;k++)for(int o=(1<<8)-1^k;~o;o=o?o-1&((1<<8)-1^k):-1)
(dp1[j+1&1][k][o]+=dp1[j&1][k][o])%=mod,(dp1[j+1&1][k|hav[i][j]][o]+=dp1[j&1][k][o])%=mod,
(dp2[j+1&1][k][o]+=dp2[j&1][k][o])%=mod,(dp2[j+1&1][k][o|hav[i][j]]+=dp2[j&1][k][o])%=mod;
}
for(int j=0;j<1<<8;j++)for(int k=(1<<8)-1^j;~k;k=k?k-1&((1<<8)-1^j):-1)
dp[now&1][j][k]=((dp1[hav[i].size()&1][j][k]+dp2[hav[i].size()&1][j][k]-dp[now-1&1][j][k])%mod+mod)%mod;
}
int ans=0;
for(int i=0;i<1<<8;i++)for(int j=(1<<8)-1^i;~j;j=j?j-1&((1<<8)-1^i):-1)(ans+=dp[now&1][i][j])%=mod;
cout<<ans;
return 0;
}
