我记得这是我 NOI 的某天晚上去自习室写的题，然后题解看不懂然后没电了，回寝室之后 AC 了，AC 了然后还是不会（
还有，洛谷咋这么喜欢关题解入口啊？？？
洛谷题目页面传送门
给定 (n)，求有多少个集合有序二元组 ((S,T)) 满足 (S,Tsubseteq{2,3,cdots,n}) 且 (forall xin S,forall yin T,xperp y)。答案对 (p) 取模。
(nin[2,500],pinleft[1,10^{10} ight])。
显然，不能打表
显然，条件可以转化为：(S,T) 分别的质因数集合没有交集。于是很容易想到状压 DP。
容易想到 (2) 个比较暴力的状压 DP：
- (dp_{i,j}) 表示考虑到数 (i)，(S) 的质因数集合为 (j) 的方案数。很好转移，刷表即可；比较难的地方在于最后统计答案，注意到 (S,T) 的情况应该是一样的，所以答案就是 (sumlimits_{Acap B=varnothing}dp_{n,A}dp_{n,B})。时间复杂度 (mathrm O!left(n2^{pi(n)} ight))。
- 不像上面那样绕弯子了，直白一点：(dp_{i,j,k}) 表示考虑到数 (i)，(S,T) 的质因数集合分别为 (j,k) 的方案数。转移依然刷表；统计答案就直接统计了。时间复杂度 (mathrm O!left(n4^{pi(n)} ight))。
写个程序随便算一下发现 (n=500) 时 (pi(n)=95)，无论哪种都是跑不过的。
注意到一个性质：一个数的所有质因数之积要小于等于原数。那么显然，(geqsqrt 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)}
ight)) 的 DP 数组，毛估估一下理论可过，但是 (mod) 运算写的丑的话会被卡常，其实可以优化。空间可以滚动数组滚成总 (mathrm O!left(4^{pi(sqrt n)}
ight)) 的；时间上的话，注意到若 DP 数组两个维度有交集的话，一定是不合法的，直接不看，可以优化到 (mathrm O!left(n3^{pi(sqrt n)}
ight))，但我比较懒，初始化的时候直接 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;
}