好题呀。
首先两者所拿的集合必然是不相交的。首先我们考虑30分的暴力做法，30以内质数很少，我们可以进行一波状压，用(dp[i][j])表示第一个人取集合i，第二个人取集合j的方案数。我们记录当前数字的质因子集合为k，之后就可以进行转移了。
但是数据范围到了500，就会导致无法状压那么多的质数，我们考虑500以内的数，最多只可能有1个大于22的质因子。
这就意味着我们完全可以把它单独提出来处理，之后剩下的8个质数我们继续进行状压。
具体的操作方法就是，我们首先把所有的数按照最大质因子进行排序，相同的排在一起，因为相同的两者是不能同时选的。对于每个最大质因子，我们用(f1[][])表示这个最大质因子让第一个人选取，(f2[][])表示这个最大质因子让第二个人选取。
之后在小范围内进行状压DP即可。
记录当前元素的质因子集合为k，转移方程是:
(f1[S_1 | k][S_2] += f1[S_1][S_2])((k & S_2 == 0))
(f2[S_1][S_2 | k] += f2[S_1][S_2])((k & S_1 == 0))
之后合并答案的时候，(dp[i][j] = f1[i][j] + f2[i][j] - dp[i][j])，其中减去一次是考虑到两个人都没有选择这个最大质因子，这种情况被计算了两次。
时间复杂度(O(n2^{16}))
#include<bits/stdc++.h>
#define rep(i,a,n) for(int i = a;i <= n;i++)
#define per(i,n,a) for(int i = n;i >= a;i--)
#define enter putchar('
')
#define fr friend inline
#define y1 poj
using namespace std;
typedef long long ll;
const int M = 1005;
const int INF = 1000000009;
const double eps = 1e-7;
int read()
{
int ans = 0,op = 1;char ch = getchar();
while(ch < '0' || ch > '9') {if(ch == '-') op = -1;ch = getchar();}
while(ch >= '0' && ch <= '9') ans = ans * 10 + ch - '0',ch = getchar();
return ans * op;
}
int n,P,p[9] = {0,2,3,5,7,11,13,17,19};
int tot,dp[300][300],f1[300][300],f2[300][300],ans;
bool np[M];
struct node
{
int val,mprime,S;
bool operator < (const node &g){return mprime < g.mprime;}
}a[M];
int inc(int a,int b){return (a+b) % P;}
void get(node &x)
{
int cur = 0,t = x.val;
rep(i,1,8)
{
if(t % p[i]) continue;
cur |= (1 << (i-1));
while(!(t % p[i])) t /= p[i];
}
x.S = cur,x.mprime = (t > 1) ? t : -1;
}
int main()
{
n = read(),P = read();
rep(i,2,n) a[i-1].val = i,get(a[i-1]);
sort(a+1,a+n),dp[0][0] = 1;
rep(i,1,n-1)
{
if(i == 1 || a[i].mprime != a[i-1].mprime || a[i].mprime == -1) memcpy(f1,dp,sizeof(f1)),memcpy(f2,dp,sizeof(f2));
per(j,255,0)
per(k,255,0)
{
if(j & k) continue;
if(!(a[i].S & j)) f2[j][k | a[i].S] = inc(f2[j][k | a[i].S],f2[j][k]);
if(!(a[i].S & k)) f1[j | a[i].S][k] = inc(f1[j | a[i].S][k],f1[j][k]);
}
if(i == n-1 || a[i].mprime != a[i+1].mprime || a[i].mprime == -1)
{
rep(j,0,255)
rep(k,0,255) if(!(j & k)) dp[j][k] = inc(f1[j][k],inc(f2[j][k],P - dp[j][k]));
}
}
rep(i,0,255)
rep(j,0,255) if(!(i&j) && dp[i][j]) ans = inc(ans,dp[i][j]);
printf("%d
",ans);
return 0;
}